LCOV - code coverage report
Current view: top level - src/motion - neb_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 389 435 89.4 %
Date: 2024-04-20 06:29:22 Functions: 10 11 90.9 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Module with utility 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_utils
      21             :    USE bibliography,                    ONLY: E2002,&
      22             :                                               Elber1987,&
      23             :                                               Jonsson1998,&
      24             :                                               Jonsson2000_1,&
      25             :                                               Jonsson2000_2,&
      26             :                                               Wales2004,&
      27             :                                               cite_reference
      28             :    USE colvar_utils,                    ONLY: eval_colvar,&
      29             :                                               get_clv_force,&
      30             :                                               set_colvars_target
      31             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      32             :                                               cp_logger_type
      33             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      34             :                                               cp_print_key_unit_nr
      35             :    USE cp_parser_methods,               ONLY: parser_get_next_line,&
      36             :                                               parser_get_object
      37             :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      38             :                                               parser_create,&
      39             :                                               parser_release
      40             :    USE f77_interface,                   ONLY: f_env_add_defaults,&
      41             :                                               f_env_rm_defaults,&
      42             :                                               f_env_type,&
      43             :                                               get_energy,&
      44             :                                               get_force,&
      45             :                                               get_pos,&
      46             :                                               set_pos
      47             :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      48             :    USE force_env_types,                 ONLY: force_env_get
      49             :    USE geo_opt,                         ONLY: cp_geo_opt
      50             :    USE global_types,                    ONLY: global_environment_type
      51             :    USE input_constants,                 ONLY: &
      52             :         band_diis_opt, do_b_neb, do_band_cartesian, do_ci_neb, do_d_neb, do_eb, do_it_neb, do_sm, &
      53             :         pot_neb_fe, pot_neb_full, pot_neb_me
      54             :    USE input_cp2k_check,                ONLY: remove_restart_info
      55             :    USE input_section_types,             ONLY: section_vals_get,&
      56             :                                               section_vals_get_subs_vals,&
      57             :                                               section_vals_type,&
      58             :                                               section_vals_val_get
      59             :    USE kinds,                           ONLY: default_path_length,&
      60             :                                               default_string_length,&
      61             :                                               dp
      62             :    USE md_run,                          ONLY: qs_mol_dyn
      63             :    USE message_passing,                 ONLY: mp_para_env_type
      64             :    USE neb_io,                          ONLY: dump_replica_coordinates,&
      65             :                                               handle_band_file_names
      66             :    USE neb_md_utils,                    ONLY: neb_initialize_velocity
      67             :    USE neb_types,                       ONLY: neb_type,&
      68             :                                               neb_var_type
      69             :    USE particle_types,                  ONLY: particle_type
      70             :    USE physcon,                         ONLY: bohr
      71             :    USE replica_methods,                 ONLY: rep_env_calc_e_f
      72             :    USE replica_types,                   ONLY: rep_env_sync,&
      73             :                                               replica_env_type
      74             :    USE rmsd,                            ONLY: rmsd3
      75             : #include "../base/base_uses.f90"
      76             : 
      77             :    IMPLICIT NONE
      78             :    PRIVATE
      79             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_utils'
      80             :    LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
      81             : 
      82             :    PUBLIC :: build_replica_coords, &
      83             :              neb_calc_energy_forces, &
      84             :              reorient_images, &
      85             :              reparametrize_images, &
      86             :              check_convergence
      87             : 
      88             : CONTAINS
      89             : 
      90             : ! **************************************************************************************************
      91             : !> \brief Computes the distance between two replica
      92             : !> \param particle_set ...
      93             : !> \param coords ...
      94             : !> \param i0 ...
      95             : !> \param i ...
      96             : !> \param distance ...
      97             : !> \param iw ...
      98             : !> \param rotate ...
      99             : !> \author Teodoro Laino 09.2006
     100             : ! **************************************************************************************************
     101       13686 :    SUBROUTINE neb_replica_distance(particle_set, coords, i0, i, distance, iw, rotate)
     102             :       TYPE(particle_type), DIMENSION(:), OPTIONAL, &
     103             :          POINTER                                         :: particle_set
     104             :       TYPE(neb_var_type), POINTER                        :: coords
     105             :       INTEGER, INTENT(IN)                                :: i0, i
     106             :       REAL(KIND=dp), INTENT(OUT)                         :: distance
     107             :       INTEGER, INTENT(IN)                                :: iw
     108             :       LOGICAL, INTENT(IN), OPTIONAL                      :: rotate
     109             : 
     110             :       LOGICAL                                            :: my_rotate
     111             : 
     112       13686 :       my_rotate = .FALSE.
     113       13686 :       IF (PRESENT(rotate)) my_rotate = rotate
     114             :       ! The rotation of the replica is enabled exclusively when working in
     115             :       ! cartesian coordinates
     116       13686 :       IF (my_rotate .AND. (coords%in_use == do_band_cartesian)) THEN
     117         444 :          CPASSERT(PRESENT(particle_set))
     118             :          CALL rmsd3(particle_set, coords%xyz(:, i), coords%xyz(:, i0), &
     119         444 :                     iw, rotate=my_rotate)
     120             :       END IF
     121             :       distance = SQRT(DOT_PRODUCT(coords%wrk(:, i) - coords%wrk(:, i0), &
     122     1853878 :                                   coords%wrk(:, i) - coords%wrk(:, i0)))
     123             : 
     124       13686 :    END SUBROUTINE neb_replica_distance
     125             : 
     126             : ! **************************************************************************************************
     127             : !> \brief Constructs or Read the coordinates for all replica
     128             : !> \param neb_section ...
     129             : !> \param particle_set ...
     130             : !> \param coords ...
     131             : !> \param vels ...
     132             : !> \param neb_env ...
     133             : !> \param iw ...
     134             : !> \param globenv ...
     135             : !> \param para_env ...
     136             : !> \author Teodoro Laino 09.2006
     137             : ! **************************************************************************************************
     138          34 :    SUBROUTINE build_replica_coords(neb_section, particle_set, &
     139             :                                    coords, vels, neb_env, iw, globenv, para_env)
     140             :       TYPE(section_vals_type), POINTER                   :: neb_section
     141             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     142             :       TYPE(neb_var_type), POINTER                        :: coords, vels
     143             :       TYPE(neb_type), POINTER                            :: neb_env
     144             :       INTEGER, INTENT(IN)                                :: iw
     145             :       TYPE(global_environment_type), POINTER             :: globenv
     146             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     147             : 
     148             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_replica_coords'
     149             : 
     150             :       CHARACTER(LEN=default_path_length)                 :: filename
     151             :       INTEGER :: handle, i_rep, iatom, ic, input_nr_replica, is, ivar, j, jtarg, n_rep, natom, &
     152             :          neb_nr_replica, nr_replica_to_interpolate, nval, nvar, shell_index
     153             :       LOGICAL                                            :: check, explicit, skip_vel_section
     154          34 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: distance
     155             :       REAL(KIND=dp), DIMENSION(3)                        :: r
     156          34 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: initial_colvars, rptr
     157             :       TYPE(section_vals_type), POINTER                   :: coord_section, replica_section, &
     158             :                                                             vel_section
     159             : 
     160          34 :       CALL timeset(routineN, handle)
     161          34 :       CPASSERT(ASSOCIATED(coords))
     162          34 :       CPASSERT(ASSOCIATED(vels))
     163          34 :       neb_nr_replica = neb_env%number_of_replica
     164          34 :       replica_section => section_vals_get_subs_vals(neb_section, "REPLICA")
     165          34 :       CALL section_vals_get(replica_section, n_repetition=input_nr_replica)
     166             :       ! Calculation is aborted if input replicas are more then the requested ones for the BAND..
     167          34 :       CPASSERT(input_nr_replica <= neb_nr_replica)
     168             :       ! Read in replicas coordinates
     169          34 :       skip_vel_section = (input_nr_replica /= neb_nr_replica)
     170          34 :       IF ((iw > 0) .AND. skip_vel_section) THEN
     171           2 :          WRITE (iw, '(T2,A)') 'NEB| The number of replica in the input is different from the number', &
     172           2 :             'NEB| of replica requested for NEB. More Replica will be interpolated.', &
     173           4 :             'NEB| Therefore the possibly provided velocities will not be read.'
     174             :       END IF
     175             :       ! Further check on velocity section...
     176         194 :       DO i_rep = 1, input_nr_replica
     177             :          vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", &
     178         160 :                                                    i_rep_section=i_rep)
     179         160 :          CALL section_vals_get(vel_section, explicit=explicit)
     180         222 :          skip_vel_section = skip_vel_section .OR. (.NOT. explicit)
     181             :       END DO
     182             :       ! Setup cartesian coordinates and COLVAR (if requested)
     183       41772 :       coords%xyz(:, :) = 0.0_dp
     184         194 :       DO i_rep = 1, input_nr_replica
     185             :          coord_section => section_vals_get_subs_vals(replica_section, "COORD", &
     186         160 :                                                      i_rep_section=i_rep)
     187         160 :          CALL section_vals_get(coord_section, explicit=explicit)
     188             :          ! Cartesian Coordinates
     189         160 :          IF (explicit) THEN
     190             :             CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", &
     191          32 :                                       n_rep_val=natom)
     192          32 :             CPASSERT((natom == SIZE(particle_set)))
     193        1512 :             DO iatom = 1, natom
     194             :                CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", &
     195        1480 :                                          i_rep_val=iatom, r_vals=rptr)
     196        1480 :                ic = 3*(iatom - 1)
     197       10360 :                coords%xyz(ic + 1:ic + 3, i_rep) = rptr(1:3)*bohr
     198             :                ! Initially core and shell positions are set to the atomic positions
     199        1480 :                shell_index = particle_set(iatom)%shell_index
     200        1512 :                IF (shell_index /= 0) THEN
     201        1140 :                   is = 3*(natom + shell_index - 1)
     202        7980 :                   coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep)
     203             :                END IF
     204             :             END DO
     205             :          ELSE
     206         128 :             BLOCK
     207             :                LOGICAL :: my_end
     208             :                CHARACTER(LEN=default_string_length)               :: dummy_char
     209             :                TYPE(cp_parser_type) :: parser
     210             :                CALL section_vals_val_get(replica_section, "COORD_FILE_NAME", &
     211         128 :                                          i_rep_section=i_rep, c_val=filename)
     212         128 :                CPASSERT(TRIM(filename) /= "")
     213         128 :                CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.TRUE.)
     214         128 :                CALL parser_get_next_line(parser, 1)
     215             :                ! Start parser
     216         128 :                CALL parser_get_object(parser, natom)
     217         128 :                CPASSERT((natom == SIZE(particle_set)))
     218         128 :                CALL parser_get_next_line(parser, 1)
     219        9774 :                DO iatom = 1, natom
     220             :                   ! Atom coordinates
     221        9646 :                   CALL parser_get_next_line(parser, 1, at_end=my_end)
     222        9646 :                   IF (my_end) &
     223             :                      CALL cp_abort(__LOCATION__, &
     224             :                                    "Number of lines in XYZ format not equal to the number of atoms."// &
     225             :                                    " Error in XYZ format for REPLICA coordinates. Very probably the"// &
     226           0 :                                    " line with title is missing or is empty. Please check the XYZ file and rerun your job!")
     227        9646 :                   READ (parser%input_line, *) dummy_char, r(1:3)
     228        9646 :                   ic = 3*(iatom - 1)
     229       38584 :                   coords%xyz(ic + 1:ic + 3, i_rep) = r(1:3)*bohr
     230             :                   ! Initially core and shell positions are set to the atomic positions
     231        9646 :                   shell_index = particle_set(iatom)%shell_index
     232        9774 :                   IF (shell_index /= 0) THEN
     233           0 :                      is = 3*(natom + shell_index - 1)
     234           0 :                      coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep)
     235             :                   END IF
     236             :                END DO
     237         512 :                CALL parser_release(parser)
     238             :             END BLOCK
     239             :          END IF
     240             :          ! Collective Variables
     241         160 :          IF (neb_env%use_colvar) THEN
     242             :             CALL section_vals_val_get(replica_section, "COLLECTIVE", &
     243          18 :                                       i_rep_section=i_rep, n_rep_val=n_rep)
     244          18 :             IF (n_rep /= 0) THEN
     245             :                ! Read the values of the collective variables
     246          10 :                NULLIFY (initial_colvars)
     247             :                CALL section_vals_val_get(replica_section, "COLLECTIVE", &
     248          10 :                                          i_rep_section=i_rep, r_vals=initial_colvars)
     249          10 :                check = (neb_env%nsize_int == SIZE(initial_colvars))
     250          10 :                CPASSERT(check)
     251          30 :                coords%int(:, i_rep) = initial_colvars
     252             :             ELSE
     253             :                ! Compute the values of the collective variables
     254           8 :                CALL eval_colvar(neb_env%force_env, coords%xyz(:, i_rep), coords%int(:, i_rep))
     255             :             END IF
     256             :          END IF
     257             :          ! Dump cartesian and colvar info..
     258         160 :          CALL dump_replica_coordinates(particle_set, coords, i_rep, i_rep, iw, neb_env%use_colvar)
     259             :          ! Setup Velocities
     260         354 :          IF (skip_vel_section) THEN
     261             :             CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, &
     262         132 :                                          i_rep, iw, globenv, neb_env)
     263             :          ELSE
     264             :             vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", &
     265          28 :                                                       i_rep_section=i_rep)
     266             :             CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
     267          28 :                                       n_rep_val=nval)
     268             :             ! Setup Velocities for collective or cartesian coordinates
     269          28 :             IF (neb_env%use_colvar) THEN
     270          10 :                nvar = SIZE(vels%wrk, 1)
     271          10 :                CPASSERT(nval == nvar)
     272          20 :                DO ivar = 1, nvar
     273             :                   CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
     274          10 :                                             i_rep_val=ivar, r_vals=rptr)
     275          20 :                   vels%wrk(ivar, i_rep) = rptr(1)
     276             :                END DO
     277             :             ELSE
     278          18 :                natom = SIZE(particle_set)
     279          18 :                CPASSERT(nval == natom)
     280         948 :                DO iatom = 1, natom
     281             :                   CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
     282         930 :                                             i_rep_val=iatom, r_vals=rptr)
     283         930 :                   ic = 3*(iatom - 1)
     284        6510 :                   vels%wrk(ic + 1:ic + 3, i_rep) = rptr(1:3)
     285             :                   ! Initially set shell velocities to core velocity
     286         930 :                   shell_index = particle_set(iatom)%shell_index
     287         948 :                   IF (shell_index /= 0) THEN
     288         760 :                      is = 3*(natom + shell_index - 1)
     289        5320 :                      vels%wrk(is + 1:is + 3, i_rep) = vels%wrk(ic + 1:ic + 3, i_rep)
     290             :                   END IF
     291             :                END DO
     292             :             END IF
     293             :          END IF
     294             :       END DO ! i_rep
     295         102 :       ALLOCATE (distance(neb_nr_replica - 1))
     296          34 :       IF (input_nr_replica < neb_nr_replica) THEN
     297             :          ! Interpolate missing replicas
     298          12 :          nr_replica_to_interpolate = neb_nr_replica - input_nr_replica
     299         104 :          distance = 0.0_dp
     300          12 :          IF (iw > 0) THEN
     301           2 :             WRITE (iw, '(T2,A,I0,A)') 'NEB| Interpolating ', nr_replica_to_interpolate, ' missing Replica.'
     302             :          END IF
     303          64 :          DO WHILE (nr_replica_to_interpolate > 0)
     304             :             ! Compute distances between known images to find the interval
     305             :             ! where to add a new image
     306         358 :             DO j = 1, input_nr_replica - 1
     307             :                CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
     308         358 :                                          rotate=neb_env%align_frames)
     309             :             END DO
     310         410 :             jtarg = MAXLOC(distance(1:input_nr_replica), 1)
     311          52 :             IF (iw > 0) THEN
     312           3 :                WRITE (iw, '(T2,3(A,I0),A)') 'NEB| Interpolating Nr. ', &
     313           3 :                   nr_replica_to_interpolate, ' missing Replica between Replica Nr. ', &
     314           6 :                   jtarg, ' and ', jtarg + 1, '.'
     315             :             END IF
     316          52 :             input_nr_replica = input_nr_replica + 1
     317          52 :             nr_replica_to_interpolate = nr_replica_to_interpolate - 1
     318             :             ! Interpolation is a simple bisection in XYZ
     319       12422 :             coords%xyz(:, jtarg + 2:input_nr_replica) = coords%xyz(:, jtarg + 1:input_nr_replica - 1)
     320        4780 :             coords%xyz(:, jtarg + 1) = (coords%xyz(:, jtarg) + coords%xyz(:, jtarg + 2))/2.0_dp
     321          52 :             IF (neb_env%use_colvar) THEN
     322             :                ! Interpolation is a simple bisection also in internal coordinates
     323             :                ! in this case the XYZ coordinates need only as a starting point for computing
     324             :                ! the potential energy function. The reference are the internal coordinates
     325             :                ! interpolated here after..
     326           6 :                coords%int(:, jtarg + 2:input_nr_replica) = coords%int(:, jtarg + 1:input_nr_replica - 1)
     327           4 :                coords%int(:, jtarg + 1) = (coords%int(:, jtarg) + coords%int(:, jtarg + 2))/2.0_dp
     328             :             END IF
     329       12322 :             vels%wrk(:, jtarg + 2:input_nr_replica) = vels%wrk(:, jtarg + 1:input_nr_replica - 1)
     330        4680 :             vels%wrk(:, jtarg + 1) = 0.0_dp
     331             :             CALL dump_replica_coordinates(particle_set, coords, jtarg + 1, &
     332          52 :                                           input_nr_replica, iw, neb_env%use_colvar)
     333             :             CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, &
     334          64 :                                          jtarg + 1, iw, globenv, neb_env)
     335             :          END DO
     336             :       END IF
     337        8126 :       vels%wrk(:, 1) = 0.0_dp
     338        8126 :       vels%wrk(:, neb_nr_replica) = 0.0_dp
     339             :       ! If we perform a DIIS optimization we don't need velocities
     340       37092 :       IF (neb_env%opt_type == band_diis_opt) vels%wrk = 0.0_dp
     341             :       ! Compute distances between replicas and in case of Cartesian Coordinates
     342             :       ! Rotate the frames in order to minimize the RMSD
     343         212 :       DO j = 1, input_nr_replica - 1
     344             :          CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
     345         212 :                                    rotate=neb_env%align_frames)
     346             :       END DO
     347          34 :       DEALLOCATE (distance)
     348             : 
     349          34 :       CALL timestop(handle)
     350             : 
     351          68 :    END SUBROUTINE build_replica_coords
     352             : 
     353             : ! **************************************************************************************************
     354             : !> \brief Driver to compute energy and forces within a NEB,
     355             : !>      Based on the use of the replica_env
     356             : !> \param rep_env ...
     357             : !> \param neb_env ...
     358             : !> \param coords ...
     359             : !> \param energies ...
     360             : !> \param forces ...
     361             : !> \param particle_set ...
     362             : !> \param output_unit ...
     363             : !> \author Teodoro Laino 09.2006
     364             : ! **************************************************************************************************
     365         908 :    SUBROUTINE neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, &
     366             :                                      particle_set, output_unit)
     367             :       TYPE(replica_env_type), POINTER                    :: rep_env
     368             :       TYPE(neb_type), OPTIONAL, POINTER                  :: neb_env
     369             :       TYPE(neb_var_type), POINTER                        :: coords
     370             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: energies
     371             :       TYPE(neb_var_type), POINTER                        :: forces
     372             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     373             :       INTEGER, INTENT(IN)                                :: output_unit
     374             : 
     375             :       CHARACTER(len=*), PARAMETER :: routineN = 'neb_calc_energy_forces'
     376             :       CHARACTER(LEN=1), DIMENSION(3), PARAMETER          :: lab = (/"X", "Y", "Z"/)
     377             : 
     378             :       INTEGER                                            :: handle, i, irep, j, n_int, n_rep, &
     379             :                                                             n_rep_neb, nsize_wrk
     380         908 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: tangent, tmp_a, tmp_b
     381             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cvalues, Mmatrix, Mmatrix_tmp
     382             : 
     383         908 :       CALL timeset(routineN, handle)
     384         908 :       n_int = neb_env%nsize_int
     385         908 :       n_rep_neb = neb_env%number_of_replica
     386         908 :       n_rep = rep_env%nrep
     387         908 :       nsize_wrk = coords%size_wrk(1)
     388        5948 :       energies = 0.0_dp
     389        2748 :       ALLOCATE (cvalues(n_int, n_rep))
     390        2748 :       ALLOCATE (Mmatrix_tmp(n_int*n_int, n_rep))
     391        2748 :       ALLOCATE (Mmatrix(n_int*n_int, n_rep_neb))
     392         908 :       IF (output_unit > 0) WRITE (output_unit, '(/,T2,A)') "NEB| Computing Energies and Forces"
     393        3700 :       DO irep = 1, n_rep_neb, n_rep
     394        8346 :          DO j = 0, n_rep - 1
     395        8346 :             IF (irep + j <= n_rep_neb) THEN
     396             :                ! If the number of replica in replica_env and the number of replica
     397             :                ! used in the NEB does not match, the other replica in replica_env
     398             :                ! just compute energies and forces keeping the fixed coordinates and
     399             :                ! forces
     400     2178252 :                rep_env%r(:, j + 1) = coords%xyz(:, irep + j)
     401             :             END IF
     402             :          END DO
     403             :          ! Fix file name for BAND replicas.. Each BAND replica has its own file
     404             :          ! independently from the number of replicas in replica_env..
     405        2792 :          CALL handle_band_file_names(rep_env, irep, n_rep_neb, neb_env%istep)
     406             :          ! Let's select the potential we want to use for the band calculation
     407        5512 :          SELECT CASE (neb_env%pot_type)
     408             :          CASE (pot_neb_full)
     409             :             ! Full potential Energy
     410        2720 :             CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
     411             :          CASE (pot_neb_fe)
     412             :             ! Free Energy Case
     413           0 :             CALL perform_replica_md(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix_tmp)
     414             :          CASE (pot_neb_me)
     415             :             ! Minimum Potential Energy Case
     416        2792 :             CALL perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix_tmp)
     417             :          END SELECT
     418             : 
     419        9254 :          DO j = 0, n_rep - 1
     420        8346 :             IF (irep + j <= n_rep_neb) THEN
     421             :                ! Copy back Forces and Energies
     422     2166252 :                forces%wrk(:, irep + j) = rep_env%f(1:nsize_wrk, j + 1)
     423        5040 :                energies(irep + j) = rep_env%f(rep_env%ndim + 1, j + 1)
     424        9960 :                SELECT CASE (neb_env%pot_type)
     425             :                CASE (pot_neb_full)
     426             :                   ! Dump Info
     427        4920 :                   IF (output_unit > 0) THEN
     428             :                      WRITE (output_unit, '(T2,A,I5,A,I5,A)') &
     429           0 :                         "NEB| REPLICA Nr.", irep + j, "- Energy and Forces"
     430             :                      WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
     431           0 :                         "NEB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
     432           0 :                      WRITE (output_unit, '(T2,"NEB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
     433           0 :                      DO i = 1, SIZE(particle_set)
     434             :                         WRITE (output_unit, '(T2,"NEB|",T12,A,T30,3(2X,F15.9))') &
     435           0 :                            particle_set(i)%atomic_kind%name, &
     436           0 :                            rep_env%f((i - 1)*3 + 1:(i - 1)*3 + 3, j + 1)
     437             :                      END DO
     438             :                   END IF
     439             :                CASE (pot_neb_fe, pot_neb_me)
     440             :                   ! Let's update the cartesian coordinates. This will make
     441             :                   ! easier the next evaluation of energies and forces...
     442       12480 :                   coords%xyz(:, irep + j) = rep_env%r(1:rep_env%ndim, j + 1)
     443         240 :                   Mmatrix(:, irep + j) = Mmatrix_tmp(:, j + 1)
     444        5160 :                   IF (output_unit > 0) THEN
     445             :                      WRITE (output_unit, '(/,T2,A,I5,A,I5,A)') &
     446          60 :                         "NEB| REPLICA Nr.", irep + j, "- Energy, Collective Variables,  Forces"
     447             :                      WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
     448          60 :                         "NEB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
     449             :                      WRITE (output_unit, &
     450          60 :                             '(T2,"NEB|",T10,"CV Nr.",12X,"Expected COLVAR",5X,"Present COLVAR",10X,"Forces")')
     451         120 :                      DO i = 1, n_int
     452             :                         WRITE (output_unit, '(T2,"NEB|",T12,I2,7X,3(5X,F15.9))') &
     453         120 :                            i, coords%int(i, irep + j), cvalues(i, j + 1), rep_env%f(i, j + 1)
     454             :                      END DO
     455             :                   END IF
     456             :                END SELECT
     457             :             END IF
     458             :          END DO
     459             :       END DO
     460         908 :       DEALLOCATE (cvalues)
     461         908 :       DEALLOCATE (Mmatrix_tmp)
     462         908 :       IF (PRESENT(neb_env)) THEN
     463             :          ! First identify the image of the chain with the higher potential energy
     464             :          ! First and last point of the band are never considered
     465        4132 :          neb_env%nr_HE_image = MAXLOC(energies(2:n_rep_neb - 1), 1) + 1
     466        2724 :          ALLOCATE (tangent(nsize_wrk))
     467             :          ! Then modify image forces accordingly to the scheme chosen for the
     468             :          ! calculation.
     469         908 :          neb_env%spring_energy = 0.0_dp
     470         908 :          IF (neb_env%optimize_end_points) THEN
     471        1290 :             ALLOCATE (tmp_a(SIZE(forces%wrk, 1)))
     472        1290 :             ALLOCATE (tmp_b(SIZE(forces%wrk, 1)))
     473      212314 :             tmp_a(:) = forces%wrk(:, 1)
     474      212314 :             tmp_b(:) = forces%wrk(:, SIZE(forces%wrk, 2))
     475             :          END IF
     476        5040 :          DO i = 2, neb_env%number_of_replica
     477        4132 :             CALL get_tangent(neb_env, coords, i, tangent, energies, output_unit)
     478             :             CALL get_neb_force(neb_env, tangent, coords, i, forces, Mmatrix=Mmatrix, &
     479        5040 :                                iw=output_unit)
     480             :          END DO
     481         908 :          IF (neb_env%optimize_end_points) THEN
     482      212314 :             forces%wrk(:, 1) = tmp_a ! Image A
     483      212314 :             forces%wrk(:, SIZE(forces%wrk, 2)) = tmp_b ! Image B
     484         430 :             DEALLOCATE (tmp_a)
     485         430 :             DEALLOCATE (tmp_b)
     486             :          ELSE
     487             :             ! Nullify forces on the two end points images
     488       37102 :             forces%wrk(:, 1) = 0.0_dp ! Image A
     489       37102 :             forces%wrk(:, SIZE(forces%wrk, 2)) = 0.0_dp ! Image B
     490             :          END IF
     491         908 :          DEALLOCATE (tangent)
     492             :       END IF
     493         908 :       DEALLOCATE (Mmatrix)
     494         908 :       CALL timestop(handle)
     495        1816 :    END SUBROUTINE neb_calc_energy_forces
     496             : 
     497             : ! **************************************************************************************************
     498             : !> \brief Driver to perform an MD run on each single replica to
     499             : !>      compute specifically Free Energies in a NEB scheme
     500             : !> \param rep_env ...
     501             : !> \param coords ...
     502             : !> \param irep ...
     503             : !> \param n_rep_neb ...
     504             : !> \param cvalues ...
     505             : !> \param Mmatrix ...
     506             : !> \author Teodoro Laino  01.2007
     507             : ! **************************************************************************************************
     508           0 :    SUBROUTINE perform_replica_md(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix)
     509             :       TYPE(replica_env_type), POINTER                    :: rep_env
     510             :       TYPE(neb_var_type), POINTER                        :: coords
     511             :       INTEGER, INTENT(IN)                                :: irep, n_rep_neb
     512             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: cvalues, Mmatrix
     513             : 
     514             :       CHARACTER(len=*), PARAMETER :: routineN = 'perform_replica_md'
     515             : 
     516             :       INTEGER                                            :: handle, handle2, ierr, j, n_el
     517             :       LOGICAL                                            :: explicit
     518             :       TYPE(cp_logger_type), POINTER                      :: logger
     519             :       TYPE(f_env_type), POINTER                          :: f_env
     520             :       TYPE(global_environment_type), POINTER             :: globenv
     521             :       TYPE(section_vals_type), POINTER                   :: md_section, root_section
     522             : 
     523           0 :       CALL timeset(routineN, handle)
     524             :       CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
     525           0 :                               handle=handle2)
     526           0 :       logger => cp_get_default_logger()
     527             :       CALL force_env_get(f_env%force_env, globenv=globenv, &
     528           0 :                          root_section=root_section)
     529           0 :       j = rep_env%local_rep_indices(1) - 1
     530           0 :       n_el = 3*rep_env%nparticle
     531           0 :       Mmatrix = 0.0_dp
     532             :       ! Syncronize position on the replica procs
     533           0 :       CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
     534           0 :       CPASSERT(ierr == 0)
     535             :       !
     536           0 :       IF (irep + j <= n_rep_neb) THEN
     537           0 :          logger%iter_info%iteration(2) = irep + j
     538           0 :          CALL remove_restart_info(root_section)
     539           0 :          md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
     540           0 :          CALL section_vals_get(md_section, explicit=explicit)
     541           0 :          CPASSERT(explicit)
     542             :          ! Let's syncronize the target of Collective Variables for this run
     543           0 :          CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env)
     544             :          ! Do a molecular dynamics and get back the derivative
     545             :          ! of the free energy w.r.t. the colvar and the metric tensor
     546           0 :          CALL qs_mol_dyn(f_env%force_env, globenv=globenv)
     547             :          ! Collect the equilibrated coordinates
     548           0 :          CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
     549           0 :          CPASSERT(ierr == 0)
     550             :          ! Write he gradients in the colvar coordinates into the replica_env array
     551             :          ! and copy back also the metric tensor..
     552             :          ! work in progress..
     553           0 :          CPABORT("")
     554           0 :          rep_env%f(:, j + 1) = 0.0_dp
     555           0 :          Mmatrix = 0.0_dp
     556             :       ELSE
     557           0 :          rep_env%r(:, j + 1) = 0.0_dp
     558           0 :          rep_env%f(:, j + 1) = 0.0_dp
     559           0 :          cvalues(:, j + 1) = 0.0_dp
     560           0 :          Mmatrix(:, j + 1) = 0.0_dp
     561             :       END IF
     562           0 :       CALL rep_env_sync(rep_env, rep_env%f)
     563           0 :       CALL rep_env_sync(rep_env, rep_env%r)
     564           0 :       CALL rep_env_sync(rep_env, cvalues)
     565           0 :       CALL rep_env_sync(rep_env, Mmatrix)
     566           0 :       CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
     567           0 :       CPASSERT(ierr == 0)
     568           0 :       CALL timestop(handle)
     569           0 :    END SUBROUTINE perform_replica_md
     570             : 
     571             : ! **************************************************************************************************
     572             : !> \brief Driver to perform a GEO_OPT run on each single replica to
     573             : !>        compute specifically minimum energies in a collective variable
     574             : !>        NEB scheme
     575             : !> \param rep_env ...
     576             : !> \param coords ...
     577             : !> \param irep ...
     578             : !> \param n_rep_neb ...
     579             : !> \param cvalues ...
     580             : !> \param Mmatrix ...
     581             : !> \author Teodoro Laino 05.2007
     582             : ! **************************************************************************************************
     583          72 :    SUBROUTINE perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix)
     584             :       TYPE(replica_env_type), POINTER                    :: rep_env
     585             :       TYPE(neb_var_type), POINTER                        :: coords
     586             :       INTEGER, INTENT(IN)                                :: irep, n_rep_neb
     587             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: cvalues, Mmatrix
     588             : 
     589             :       CHARACTER(len=*), PARAMETER :: routineN = 'perform_replica_geo'
     590             : 
     591             :       INTEGER                                            :: handle, handle2, ierr, j, n_el
     592             :       LOGICAL                                            :: explicit
     593             :       TYPE(cp_logger_type), POINTER                      :: logger
     594             :       TYPE(f_env_type), POINTER                          :: f_env
     595             :       TYPE(global_environment_type), POINTER             :: globenv
     596             :       TYPE(section_vals_type), POINTER                   :: geoopt_section, root_section
     597             : 
     598          72 :       CALL timeset(routineN, handle)
     599             :       CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
     600          72 :                               handle=handle2)
     601          72 :       logger => cp_get_default_logger()
     602             :       CALL force_env_get(f_env%force_env, globenv=globenv, &
     603          72 :                          root_section=root_section)
     604          72 :       j = rep_env%local_rep_indices(1) - 1
     605          72 :       n_el = 3*rep_env%nparticle
     606         360 :       Mmatrix = 0.0_dp
     607             :       ! Syncronize position on the replica procs
     608          72 :       CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
     609          72 :       CPASSERT(ierr == 0)
     610          72 :       IF (irep + j <= n_rep_neb) THEN
     611          60 :          logger%iter_info%iteration(2) = irep + j
     612          60 :          CALL remove_restart_info(root_section)
     613          60 :          geoopt_section => section_vals_get_subs_vals(root_section, "MOTION%GEO_OPT")
     614          60 :          CALL section_vals_get(geoopt_section, explicit=explicit)
     615          60 :          CPASSERT(explicit)
     616             :          ! Let's syncronize the target of Collective Variables for this run
     617          60 :          CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env)
     618             :          ! Do a geometry optimization..
     619          60 :          CALL cp_geo_opt(f_env%force_env, globenv=globenv)
     620             :          ! Once the geometry optimization is ended let's do a single run
     621             :          ! without any constraints/restraints
     622             :          CALL force_env_calc_energy_force(f_env%force_env, &
     623          60 :                                           calc_force=.TRUE., skip_external_control=.TRUE.)
     624             :          ! Collect the optimized coordinates
     625          60 :          CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
     626          60 :          CPASSERT(ierr == 0)
     627             :          ! Collect the gradients in cartesian coordinates
     628          60 :          CALL get_force(rep_env%f_env_id, rep_env%f(1:n_el, j + 1), n_el, ierr)
     629          60 :          CPASSERT(ierr == 0)
     630             :          ! Copy the energy
     631          60 :          CALL get_energy(rep_env%f_env_id, rep_env%f(n_el + 1, j + 1), ierr)
     632          60 :          CPASSERT(ierr == 0)
     633             :          ! The gradients in the colvar coordinates
     634             :          CALL get_clv_force(f_env%force_env, rep_env%f(1:n_el, j + 1), rep_env%r(1:n_el, j + 1), &
     635          60 :                             SIZE(coords%xyz, 1), SIZE(coords%wrk, 1), cvalues(:, j + 1), Mmatrix(:, j + 1))
     636             :       ELSE
     637         624 :          rep_env%r(:, j + 1) = 0.0_dp
     638         636 :          rep_env%f(:, j + 1) = 0.0_dp
     639          24 :          cvalues(:, j + 1) = 0.0_dp
     640          24 :          Mmatrix(:, j + 1) = 0.0_dp
     641             :       END IF
     642          72 :       CALL rep_env_sync(rep_env, rep_env%f)
     643          72 :       CALL rep_env_sync(rep_env, rep_env%r)
     644          72 :       CALL rep_env_sync(rep_env, cvalues)
     645          72 :       CALL rep_env_sync(rep_env, Mmatrix)
     646          72 :       CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
     647          72 :       CPASSERT(ierr == 0)
     648          72 :       CALL timestop(handle)
     649          72 :    END SUBROUTINE perform_replica_geo
     650             : 
     651             : ! **************************************************************************************************
     652             : !> \brief Computes the tangent for point i of the NEB chain
     653             : !> \param neb_env ...
     654             : !> \param coords ...
     655             : !> \param i ...
     656             : !> \param tangent ...
     657             : !> \param energies ...
     658             : !> \param iw ...
     659             : !> \author Teodoro Laino 09.2006
     660             : ! **************************************************************************************************
     661        4132 :    SUBROUTINE get_tangent(neb_env, coords, i, tangent, energies, iw)
     662             :       TYPE(neb_type), POINTER                            :: neb_env
     663             :       TYPE(neb_var_type), POINTER                        :: coords
     664             :       INTEGER, INTENT(IN)                                :: i
     665             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: tangent
     666             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: energies
     667             :       INTEGER, INTENT(IN)                                :: iw
     668             : 
     669             :       REAL(KINd=dp)                                      :: distance0, distance1, distance2, DVmax, &
     670             :                                                             Dvmin
     671             : 
     672        4132 :       CPASSERT(ASSOCIATED(coords))
     673      833710 :       tangent(:) = 0.0_dp
     674             :       ! For the last point we don't need any tangent..
     675        4132 :       IF (i == neb_env%number_of_replica) RETURN
     676             :       ! Several kind of tangents implemented...
     677        3224 :       SELECT CASE (neb_env%id_type)
     678             :       CASE (do_eb)
     679       38792 :          tangent(:) = 0.0_dp
     680             :       CASE (do_b_neb)
     681             :          CALL neb_replica_distance(coords=coords, i0=i, i=i - 1, distance=distance1, iw=iw, &
     682         126 :                                    rotate=.FALSE.)
     683             :          CALL neb_replica_distance(coords=coords, i0=i + 1, i=i, distance=distance2, iw=iw, &
     684         126 :                                    rotate=.FALSE.)
     685             :          tangent(:) = (coords%wrk(:, i) - coords%wrk(:, i - 1))/distance1 + &
     686        6552 :                       (coords%wrk(:, i + 1) - coords%wrk(:, i))/distance2
     687             :       CASE (do_it_neb, do_ci_neb, do_d_neb)
     688        1974 :          IF ((energies(i + 1) .GT. energies(i)) .AND. (energies(i) .GT. (energies(i - 1)))) THEN
     689      226914 :             tangent(:) = coords%wrk(:, i + 1) - coords%wrk(:, i)
     690        1282 :          ELSE IF ((energies(i + 1) .LT. energies(i)) .AND. (energies(i) .LT. (energies(i - 1)))) THEN
     691       35162 :             tangent(:) = coords%wrk(:, i) - coords%wrk(:, i - 1)
     692             :          ELSE
     693        1188 :             DVmax = MAX(ABS(energies(i + 1) - energies(i)), ABS(energies(i - 1) - energies(i)))
     694        1188 :             DVmin = MIN(ABS(energies(i + 1) - energies(i)), ABS(energies(i - 1) - energies(i)))
     695        1188 :             IF (energies(i + 1) .GE. energies(i - 1)) THEN
     696       26028 :                tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*DVmax + (coords%wrk(:, i) - coords%wrk(:, i - 1))*DVmin
     697             :             ELSE
     698      231190 :                tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*DVmin + (coords%wrk(:, i) - coords%wrk(:, i - 1))*DVmax
     699             :             END IF
     700             :          END IF
     701             :       CASE (do_sm)
     702             :          ! String method..
     703       22502 :          tangent(:) = 0.0_dp
     704             :       END SELECT
     705      584294 :       distance0 = SQRT(DOT_PRODUCT(tangent(:), tangent(:)))
     706      526970 :       IF (distance0 /= 0.0_dp) tangent(:) = tangent(:)/distance0
     707             :    END SUBROUTINE get_tangent
     708             : 
     709             : ! **************************************************************************************************
     710             : !> \brief Computes the forces for point i of the NEB chain
     711             : !> \param neb_env ...
     712             : !> \param tangent ...
     713             : !> \param coords ...
     714             : !> \param i ...
     715             : !> \param forces ...
     716             : !> \param tag ...
     717             : !> \param Mmatrix ...
     718             : !> \param iw ...
     719             : !> \author Teodoro Laino 09.2006
     720             : ! **************************************************************************************************
     721        4630 :    RECURSIVE SUBROUTINE get_neb_force(neb_env, tangent, coords, i, forces, tag, Mmatrix, &
     722             :                                       iw)
     723             :       TYPE(neb_type), POINTER                            :: neb_env
     724             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: tangent
     725             :       TYPE(neb_var_type), POINTER                        :: coords
     726             :       INTEGER, INTENT(IN)                                :: i
     727             :       TYPE(neb_var_type), POINTER                        :: forces
     728             :       INTEGER, INTENT(IN), OPTIONAL                      :: tag
     729             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Mmatrix
     730             :       INTEGER, INTENT(IN)                                :: iw
     731             : 
     732             :       INTEGER                                            :: j, my_tag, nsize_wrk
     733             :       REAL(KIND=dp)                                      :: distance1, distance2, fac, tmp
     734        4630 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: dtmp1, wrk
     735             : 
     736        4630 :       my_tag = neb_env%id_type
     737        4630 :       IF (PRESENT(tag)) my_tag = tag
     738        4630 :       CPASSERT(ASSOCIATED(forces))
     739        4630 :       CPASSERT(ASSOCIATED(coords))
     740        4630 :       nsize_wrk = coords%size_wrk(1)
     741             :       ! All methods but not the classical elastic band will skip the force
     742             :       ! calculation for the last frame of the band
     743        3336 :       SELECT CASE (my_tag)
     744             :       CASE (do_b_neb, do_it_neb, do_ci_neb, do_d_neb)
     745        4494 :          IF (i == neb_env%number_of_replica) RETURN
     746             :       CASE (do_sm)
     747             :          ! String Method
     748             :          ! The forces do not require any projection. Reparametrization required
     749             :          ! after the update of the replica.
     750         420 :          CALL cite_reference(E2002)
     751        4630 :          RETURN
     752             :       END SELECT
     753             :       ! otherwise proceeed normally..
     754       10416 :       ALLOCATE (wrk(nsize_wrk))
     755             :       ! Spring Energy
     756             :       CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance1, iw=iw, &
     757        3472 :                                 rotate=.FALSE.)
     758        3472 :       tmp = distance1 - neb_env%avg_distance
     759        3472 :       neb_env%spring_energy = neb_env%spring_energy + 0.5_dp*neb_env%k*tmp**2
     760         874 :       SELECT CASE (my_tag)
     761             :       CASE (do_eb)
     762         874 :          CALL cite_reference(Elber1987)
     763             :          ! Elastic band - Hamiltonian formulation according the original Karplus/Elber
     764             :          !                formulation
     765        2622 :          ALLOCATE (dtmp1(nsize_wrk))
     766             :          ! derivatives of the spring
     767         874 :          tmp = distance1 - neb_env%avg_distance
     768       45448 :          dtmp1(:) = 1.0_dp/distance1*(coords%wrk(:, i) - coords%wrk(:, i - 1))
     769       45448 :          wrk(:) = neb_env%k*tmp*dtmp1
     770       45448 :          forces%wrk(:, i) = forces%wrk(:, i) - wrk
     771       45448 :          forces%wrk(:, i - 1) = forces%wrk(:, i - 1) + wrk
     772             :          ! derivatives due to the average length of the spring
     773         874 :          fac = 1.0_dp/(neb_env%avg_distance*REAL(neb_env%number_of_replica - 1, KIND=dp))
     774       45448 :          wrk(:) = neb_env%k*fac*(coords%wrk(:, i) - coords%wrk(:, i - 1))
     775         874 :          tmp = 0.0_dp
     776        7880 :          DO j = 2, neb_env%number_of_replica
     777             :             CALL neb_replica_distance(coords=coords, i0=j - 1, i=j, distance=distance1, iw=iw, &
     778        7006 :                                       rotate=.FALSE.)
     779        7880 :             tmp = tmp + distance1 - neb_env%avg_distance
     780             :          END DO
     781       45448 :          forces%wrk(:, i) = forces%wrk(:, i) + wrk*tmp
     782       45448 :          forces%wrk(:, i - 1) = forces%wrk(:, i - 1) - wrk*tmp
     783         874 :          DEALLOCATE (dtmp1)
     784             :       CASE (do_b_neb)
     785             :          ! Bisection NEB
     786         126 :          CALL cite_reference(Jonsson1998)
     787        6552 :          wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
     788        6552 :          tmp = neb_env%k*DOT_PRODUCT(wrk, tangent)
     789        6552 :          wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent
     790        6552 :          forces%wrk(:, i) = wrk + tmp*tangent
     791             :       CASE (do_it_neb)
     792             :          ! Improved tangent NEB
     793        1236 :          CALL cite_reference(Jonsson2000_1)
     794             :          CALL neb_replica_distance(coords=coords, i0=i, i=i + 1, distance=distance1, iw=iw, &
     795        1236 :                                    rotate=.FALSE.)
     796             :          CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance2, iw=iw, &
     797        1236 :                                    rotate=.FALSE.)
     798        1236 :          tmp = neb_env%k*(distance1 - distance2)
     799      309648 :          wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent
     800      309648 :          forces%wrk(:, i) = wrk + tmp*tangent
     801             :       CASE (do_ci_neb)
     802             :          ! Climbing Image NEB
     803         858 :          CALL cite_reference(Jonsson2000_2)
     804         858 :          IF (neb_env%istep <= neb_env%nsteps_it .OR. i /= neb_env%nr_HE_image) THEN
     805         498 :             CALL get_neb_force(neb_env, tangent, coords, i, forces, do_it_neb, Mmatrix, iw)
     806             :          ELSE
     807      189990 :             wrk(:) = forces%wrk(:, i)
     808         360 :             tmp = -2.0_dp*dot_product_band(neb_env, wrk, tangent, Mmatrix)
     809      189990 :             forces%wrk(:, i) = wrk + tmp*tangent
     810             :          END IF
     811             :       CASE (do_d_neb)
     812             :          ! Doubly NEB
     813         378 :          CALL cite_reference(Wales2004)
     814        1134 :          ALLOCATE (dtmp1(nsize_wrk))
     815       19656 :          dtmp1(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent
     816       19656 :          forces%wrk(:, i) = dtmp1
     817       19656 :          tmp = SQRT(DOT_PRODUCT(dtmp1, dtmp1))
     818       19656 :          dtmp1(:) = dtmp1(:)/tmp
     819             :          ! Project out only the spring component interfering with the
     820             :          ! orthogonal gradient of the band
     821       19656 :          wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
     822       19656 :          tmp = DOT_PRODUCT(wrk, dtmp1)
     823       19656 :          dtmp1(:) = neb_env%k*(wrk(:) - tmp*dtmp1(:))
     824       19656 :          forces%wrk(:, i) = forces%wrk(:, i) + dtmp1(:)
     825        3850 :          DEALLOCATE (dtmp1)
     826             :       END SELECT
     827        3472 :       DEALLOCATE (wrk)
     828        1158 :    END SUBROUTINE get_neb_force
     829             : 
     830             : ! **************************************************************************************************
     831             : !> \brief  Handles the dot_product when using colvar.. in this case
     832             : !>         the scalar product needs to take into account the metric
     833             : !>         tensor
     834             : !> \param neb_env ...
     835             : !> \param array1 ...
     836             : !> \param array2 ...
     837             : !> \param array3 ...
     838             : !> \return ...
     839             : !> \author Teodoro Laino 09.2006
     840             : ! **************************************************************************************************
     841        2100 :    FUNCTION dot_product_band(neb_env, array1, array2, array3) RESULT(value)
     842             :       TYPE(neb_type), POINTER                            :: neb_env
     843             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: array1, array2
     844             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: array3
     845             :       REAL(KIND=dp)                                      :: value
     846             : 
     847             :       INTEGER                                            :: nsize_int
     848             :       LOGICAL                                            :: check
     849             : 
     850        2100 :       IF (neb_env%use_colvar) THEN
     851          72 :          nsize_int = neb_env%nsize_int
     852             :          check = ((SIZE(array1) /= SIZE(array2)) .OR. &
     853             :                   (SIZE(array1) /= nsize_int) .OR. &
     854          72 :                   (SIZE(array3) /= nsize_int*nsize_int))
     855             :          ! This condition should always be satisfied..
     856           0 :          CPASSERT(check)
     857         648 :          value = DOT_PRODUCT(MATMUL(RESHAPE(array3, (/nsize_int, nsize_int/)), array1), array2)
     858             :       ELSE
     859      525702 :          value = DOT_PRODUCT(array1, array2)
     860             :       END IF
     861        2100 :    END FUNCTION dot_product_band
     862             : 
     863             : ! **************************************************************************************************
     864             : !> \brief Reorient iteratively all images of the NEB chain in order to
     865             : !>      have always the smaller RMSD between two following images
     866             : !> \param rotate_frames ...
     867             : !> \param particle_set ...
     868             : !> \param coords ...
     869             : !> \param vels ...
     870             : !> \param iw ...
     871             : !> \param distances ...
     872             : !> \param number_of_replica ...
     873             : !> \author Teodoro Laino 09.2006
     874             : ! **************************************************************************************************
     875         908 :    SUBROUTINE reorient_images(rotate_frames, particle_set, coords, vels, iw, &
     876         908 :                               distances, number_of_replica)
     877             :       LOGICAL, INTENT(IN)                                :: rotate_frames
     878             :       TYPE(particle_type), DIMENSION(:), OPTIONAL, &
     879             :          POINTER                                         :: particle_set
     880             :       TYPE(neb_var_type), POINTER                        :: coords, vels
     881             :       INTEGER, INTENT(IN)                                :: iw
     882             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: distances
     883             :       INTEGER, INTENT(IN)                                :: number_of_replica
     884             : 
     885             :       INTEGER                                            :: i, k, kind
     886             :       LOGICAL                                            :: check
     887             :       REAL(KIND=dp)                                      :: xtmp
     888             :       REAL(KIND=dp), DIMENSION(3)                        :: tmp
     889             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rot
     890             : 
     891         908 :       rot = 0.0_dp
     892         908 :       rot(1, 1) = 1.0_dp
     893         908 :       rot(2, 2) = 1.0_dp
     894         908 :       rot(3, 3) = 1.0_dp
     895        5040 :       DO i = 2, number_of_replica
     896             :          ! The rotation of the replica is enabled exclusively when working in
     897             :          ! cartesian coordinates
     898        4132 :          IF (rotate_frames .AND. (coords%in_use == do_band_cartesian)) THEN
     899             :             CALL rmsd3(particle_set, coords%xyz(:, i), coords%xyz(:, i - 1), iw, &
     900         540 :                        rotate=.TRUE., rot=rot)
     901             :             ! Rotate velocities
     902       11796 :             DO k = 1, SIZE(vels%xyz, 1)/3
     903       11256 :                kind = (k - 1)*3
     904       45024 :                tmp = vels%xyz(kind + 1:kind + 3, i)
     905       45564 :                vels%xyz(kind + 1:kind + 3, i) = MATMUL(TRANSPOSE(rot), tmp)
     906             :             END DO
     907             :          END IF
     908        5040 :          IF (PRESENT(distances)) THEN
     909        4132 :             check = SIZE(distances) == (number_of_replica - 1)
     910        4132 :             CPASSERT(check)
     911             :             xtmp = DOT_PRODUCT(coords%wrk(:, i) - coords%wrk(:, i - 1), &
     912      833710 :                                coords%wrk(:, i) - coords%wrk(:, i - 1))
     913        4132 :             distances(i - 1) = SQRT(xtmp)
     914             :          END IF
     915             :       END DO
     916         908 :    END SUBROUTINE reorient_images
     917             : 
     918             : ! **************************************************************************************************
     919             : !> \brief Reparametrization of the replica for String Method with splines
     920             : !> \param reparametrize_frames ...
     921             : !> \param spline_order ...
     922             : !> \param smoothing ...
     923             : !> \param coords ...
     924             : !> \param sline ...
     925             : !> \param distances ...
     926             : !> \author Teodoro Laino - Rodolphe Vuilleumier 09.2008
     927             : ! **************************************************************************************************
     928         402 :    SUBROUTINE reparametrize_images(reparametrize_frames, spline_order, smoothing, &
     929         402 :                                    coords, sline, distances)
     930             : 
     931             :       LOGICAL, INTENT(IN)                                :: reparametrize_frames
     932             :       INTEGER, INTENT(IN)                                :: spline_order
     933             :       REAL(KIND=dp), INTENT(IN)                          :: smoothing
     934             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: coords, sline
     935             :       REAL(KIND=dp), DIMENSION(:)                        :: distances
     936             : 
     937             :       INTEGER                                            :: i, j
     938             :       REAL(KIND=dp)                                      :: avg_distance, xtmp
     939         402 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: tmp_coords
     940             : 
     941         402 :       IF (reparametrize_frames) THEN
     942         168 :          ALLOCATE (tmp_coords(SIZE(coords, 1), SIZE(coords, 2)))
     943       24066 :          tmp_coords(:, :) = coords
     944             :          ! Smoothing
     945         420 :          DO i = 2, SIZE(coords, 2) - 1
     946             :             coords(:, i) = tmp_coords(:, i)*(1.0_dp - 2.0_dp*smoothing) + &
     947       19698 :                            tmp_coords(:, i - 1)*smoothing + tmp_coords(:, i + 1)*smoothing
     948             :          END DO
     949       48090 :          sline = coords - tmp_coords + sline
     950       24066 :          tmp_coords(:, :) = coords
     951             :          ! Reparametrization
     952          84 :          SELECT CASE (spline_order)
     953             :          CASE (1)
     954             :             ! Compute distances
     955         462 :             DO i = 2, SIZE(coords, 2)
     956       21840 :                xtmp = DOT_PRODUCT(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
     957         462 :                distances(i - 1) = SQRT(xtmp)
     958             :             END DO
     959         462 :             avg_distance = SUM(distances)/REAL(SIZE(coords, 2) - 1, KIND=dp)
     960             :             ! Redistribute frames
     961         420 :             DO i = 2, SIZE(coords, 2) - 1
     962         378 :                xtmp = 0.0_dp
     963        2022 :                DO j = 1, SIZE(coords, 2) - 1
     964        1980 :                   xtmp = xtmp + distances(j)
     965        1980 :                   IF (xtmp > avg_distance*REAL(i - 1, KIND=dp)) THEN
     966         378 :                      xtmp = (xtmp - avg_distance*REAL(i - 1, KIND=dp))/distances(j)
     967       19656 :                      coords(:, i) = (1.0_dp - xtmp)*tmp_coords(:, j + 1) + xtmp*tmp_coords(:, j)
     968             :                      EXIT
     969             :                   END IF
     970             :                END DO
     971             :             END DO
     972             :             ! Re-compute distances
     973         462 :             DO i = 2, SIZE(coords, 2)
     974       21840 :                xtmp = DOT_PRODUCT(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
     975         462 :                distances(i - 1) = SQRT(xtmp)
     976             :             END DO
     977             :          CASE DEFAULT
     978          42 :             CPWARN("String Method: Spline order greater than 1 not implemented.")
     979             :          END SELECT
     980       48090 :          sline = coords - tmp_coords + sline
     981          42 :          DEALLOCATE (tmp_coords)
     982             :       END IF
     983         402 :    END SUBROUTINE reparametrize_images
     984             : 
     985             : ! **************************************************************************************************
     986             : !> \brief Checks for convergence criteria during a NEB run
     987             : !> \param neb_env ...
     988             : !> \param Dcoords ...
     989             : !> \param forces ...
     990             : !> \return ...
     991             : !> \author Teodoro Laino 10.2006
     992             : ! **************************************************************************************************
     993        2720 :    FUNCTION check_convergence(neb_env, Dcoords, forces) RESULT(converged)
     994             :       TYPE(neb_type), POINTER                            :: neb_env
     995             :       TYPE(neb_var_type), POINTER                        :: Dcoords, forces
     996             :       LOGICAL                                            :: converged
     997             : 
     998             :       CHARACTER(LEN=3), DIMENSION(4)                     :: labels
     999             :       INTEGER                                            :: iw
    1000             :       REAL(KIND=dp)                                      :: max_dr, max_force, my_max_dr, &
    1001             :                                                             my_max_force, my_rms_dr, my_rms_force, &
    1002             :                                                             rms_dr, rms_force
    1003             :       TYPE(cp_logger_type), POINTER                      :: logger
    1004             :       TYPE(section_vals_type), POINTER                   :: cc_section
    1005             : 
    1006         544 :       NULLIFY (logger, cc_section)
    1007         544 :       logger => cp_get_default_logger()
    1008         544 :       cc_section => section_vals_get_subs_vals(neb_env%neb_section, "CONVERGENCE_CONTROL")
    1009         544 :       CALL section_vals_val_get(cc_section, "MAX_DR", r_val=max_dr)
    1010         544 :       CALL section_vals_val_get(cc_section, "MAX_FORCE", r_val=max_force)
    1011         544 :       CALL section_vals_val_get(cc_section, "RMS_DR", r_val=rms_dr)
    1012         544 :       CALL section_vals_val_get(cc_section, "RMS_FORCE", r_val=rms_force)
    1013         544 :       converged = .FALSE.
    1014        2720 :       labels = " NO"
    1015      562306 :       my_max_dr = MAXVAL(ABS(Dcoords%wrk))
    1016      562306 :       my_max_force = MAXVAL(ABS(forces%wrk))
    1017      562306 :       my_rms_dr = SQRT(SUM(Dcoords%wrk*Dcoords%wrk)/REAL(SIZE(Dcoords%wrk, 1)*SIZE(Dcoords%wrk, 2), KIND=dp))
    1018      562306 :       my_rms_force = SQRT(SUM(forces%wrk*forces%wrk)/REAL(SIZE(forces%wrk, 1)*SIZE(forces%wrk, 2), KIND=dp))
    1019         544 :       IF (my_max_dr < max_dr) labels(1) = "YES"
    1020         544 :       IF (my_max_force < max_force) labels(2) = "YES"
    1021         544 :       IF (my_rms_dr < rms_dr) labels(3) = "YES"
    1022         544 :       IF (my_rms_force < rms_force) labels(4) = "YES"
    1023         586 :       IF (ALL(labels == "YES")) converged = .TRUE.
    1024             : 
    1025             :       iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "CONVERGENCE_INFO", &
    1026         544 :                                 extension=".nebLog")
    1027         544 :       IF (iw > 0) THEN
    1028             :          ! Print convergence info
    1029         543 :          WRITE (iw, FMT='(A,A)') ' **************************************', &
    1030        1086 :             '*****************************************'
    1031             :          WRITE (iw, FMT='(1X,A,2X,F8.5,5X,"[",F8.5,"]",1X,T76,"(",A,")")') &
    1032         543 :             'RMS DISPLACEMENT =', my_rms_dr, rms_dr, labels(3), &
    1033         543 :             'MAX DISPLACEMENT =', my_max_dr, max_dr, labels(1), &
    1034         543 :             'RMS FORCE        =', my_rms_force, rms_force, labels(4), &
    1035        1086 :             'MAX FORCE        =', my_max_force, max_force, labels(2)
    1036         543 :          WRITE (iw, FMT='(A,A)') ' **************************************', &
    1037        1086 :             '*****************************************'
    1038             :       END IF
    1039             :       CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, &
    1040         544 :                                         "CONVERGENCE_INFO")
    1041         544 :    END FUNCTION check_convergence
    1042             : 
    1043          72 : END MODULE neb_utils

Generated by: LCOV version 1.15