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 CALL timestop(handle)
567 INTEGER,
INTENT(IN) :: i_rep, ienum, iw
568 LOGICAL,
INTENT(IN) :: use_colvar
571 REAL(kind=
dp),
DIMENSION(3) :: r
574 WRITE (iw,
'(/,T2,"NEB|",75("*"))')
575 WRITE (iw,
'(T2,"NEB|",1X,A,I0,A)') &
576 "Geometry for Replica Nr. ", ienum,
" in Angstrom"
577 DO iatom = 1,
SIZE(particle_set)
579 WRITE (iw,
'(T2,"NEB|",1X,A10,5X,3F15.9)') &
580 trim(particle_set(iatom)%atomic_kind%name), r(1:3)*
angstrom
583 WRITE (iw,
'(/,T2,"NEB|",1X,A10)')
"COLLECTIVE VARIABLES:"
584 WRITE (iw,
'(T2,"NEB|",16X,3F15.9)') &
585 (coords%int(j, i_rep), j=1,
SIZE(coords%int(:, :), 1))
587 WRITE (iw,
'(T2,"NEB|",75("*"))')
603 INTEGER,
INTENT(IN) :: irep, n_rep, istep
605 CHARACTER(len=*),
PARAMETER :: routinen =
'handle_band_file_names'
607 CHARACTER(LEN=default_path_length) :: output_file_path, replica_proj_name
608 INTEGER :: handle, handle2, i, ierr, j, lp, unit_nr
613 CALL timeset(routinen, handle)
617 CALL force_env_get(f_env%force_env, root_section=root_section)
618 j = irep + (rep_env%local_rep_indices(1) - 1)
620 replica_proj_name = get_replica_project_name(rep_env, n_rep, j)
621 lp = len_trim(replica_proj_name)
623 c_val=trim(replica_proj_name))
624 logger%iter_info%project_name = trim(replica_proj_name)
627 output_file_path = replica_proj_name(1:lp)//
".out"
629 c_val=trim(output_file_path))
630 IF (logger%default_global_unit_nr > 0)
THEN
631 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 skip_get_unit_number=.true.)
636 WRITE (unit=logger%default_global_unit_nr, fmt=
"(/,(T2,A79))") &
637 "*******************************************************************************", &
638 "** BAND EVALUATION OF ENERGIES AND FORCES **", &
639 "*******************************************************************************"
640 WRITE (unit=logger%default_global_unit_nr, fmt=
"(T2,A,T79,A)")
"**",
"**"
641 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 "** 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 "** Band Step Nr. :", istep,
"**"
646 WRITE (unit=logger%default_global_unit_nr, fmt=
"(T2,A79)") &
647 "*******************************************************************************"
651 SELECT CASE (f_env%force_env%in_use)
653 DO i = 1, f_env%force_env%mixed_env%ngroups
654 IF (
modulo(i - 1, f_env%force_env%mixed_env%ngroups) == &
655 f_env%force_env%mixed_env%group_distribution(f_env%force_env%mixed_env%para_env%mepos))
THEN
656 sub_logger => f_env%force_env%mixed_env%sub_logger(i)%p
657 sub_logger%iter_info%project_name = replica_proj_name(1:lp)//
"-r-"//trim(adjustl(
cp_to_string(i)))
659 unit_nr = sub_logger%default_global_unit_nr
660 IF (unit_nr > 0)
THEN
663 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 unit_number=unit_nr, skip_get_unit_number=.true.)
674 CALL timestop(handle)
686 FUNCTION get_replica_project_name(rep_env, n_rep, j)
RESULT(replica_proj_name)
688 INTEGER,
INTENT(IN) :: n_rep, j
689 CHARACTER(LEN=default_path_length) :: replica_proj_name
691 CHARACTER(LEN=default_string_length) :: padding
692 INTEGER :: i, lp, ndigits
696 replica_proj_name = rep_env%original_project_name
698 ndigits = ceiling(log10(real(n_rep + 1, kind=
dp))) - &
699 ceiling(log10(real(j + 1, kind=
dp)))
704 lp = len_trim(replica_proj_name)
705 replica_proj_name(lp + 1:len(replica_proj_name)) =
"-BAND"// &
707 END FUNCTION get_replica_project_name
721 CHARACTER(LEN=default_path_length) :: replica_proj_name
722 INTEGER :: handle2, ierr, irep, n_rep, n_rep_neb, &
727 n_rep_neb = neb_env%number_of_replica
732 output_unit = logger%default_global_unit_nr
733 IF (output_unit > 0)
THEN
734 WRITE (unit=output_unit, fmt=
'(/,(T2,A79))') &
735 "*******************************************************************************", &
736 "** MAPPING OF BAND REPLICA TO REPLICA ENV **", &
737 "*******************************************************************************"
738 WRITE (unit=output_unit, fmt=
'(T2,A,I6,T32,A,T79,A)') &
739 "** Replica Env Nr.: ", rep_env%local_rep_indices(1) - 1, &
740 "working on the following BAND replicas",
"**"
741 WRITE (unit=output_unit, fmt=
'(T2,A79)') &
744 DO irep = 1, n_rep_neb, n_rep
745 replica_proj_name = get_replica_project_name(rep_env, n_rep_neb, irep + rep_env%local_rep_indices(1) - 1)
746 IF (output_unit > 0)
THEN
747 WRITE (unit=output_unit, fmt=
'(T2,A,I6,T32,A,T79,A)') &
748 "** Band Replica Nr.: ", irep + rep_env%local_rep_indices(1) - 1, &
749 "Output available on file: "//trim(replica_proj_name)//
".out",
"**"
752 IF (output_unit > 0)
THEN
753 WRITE (unit=output_unit, fmt=
'(T2,A79)') &
755 "*******************************************************************************"
756 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