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