(git:0de0cc2)
neb_utils.F
Go to the documentation of this file.
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,&
26  wales2004,&
27  cite_reference
28  USE colvar_utils, ONLY: eval_colvar,&
32  cp_logger_type
36  parser_get_object
37  USE cp_parser_types, ONLY: cp_parser_type,&
42  f_env_type,&
43  get_energy,&
44  get_force,&
45  get_pos,&
46  set_pos
49  USE geo_opt, ONLY: cp_geo_opt
50  USE global_types, ONLY: global_environment_type
51  USE input_constants, ONLY: &
57  section_vals_type,&
59  USE kinds, ONLY: default_path_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,&
67  USE neb_types, ONLY: neb_type,&
68  neb_var_type
69  USE particle_types, ONLY: particle_type
70  USE physcon, ONLY: bohr
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, &
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  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  my_rotate = .false.
113  IF (PRESENT(rotate)) my_rotate = rotate
114  ! The rotation of the replica is enabled exclusively when working in
115  ! cartesian coordinates
116  IF (my_rotate .AND. (coords%in_use == do_band_cartesian)) THEN
117  cpassert(PRESENT(particle_set))
118  CALL rmsd3(particle_set, coords%xyz(:, i), coords%xyz(:, i0), &
119  iw, rotate=my_rotate)
120  END IF
121  distance = sqrt(dot_product(coords%wrk(:, i) - coords%wrk(:, i0), &
122  coords%wrk(:, i) - coords%wrk(:, i0)))
123 
124  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: distance
155  REAL(kind=dp), DIMENSION(3) :: r
156  REAL(kind=dp), DIMENSION(:), POINTER :: initial_colvars, rptr
157  TYPE(section_vals_type), POINTER :: coord_section, replica_section, &
158  vel_section
159 
160  CALL timeset(routinen, handle)
161  cpassert(ASSOCIATED(coords))
162  cpassert(ASSOCIATED(vels))
163  neb_nr_replica = neb_env%number_of_replica
164  replica_section => section_vals_get_subs_vals(neb_section, "REPLICA")
165  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  cpassert(input_nr_replica <= neb_nr_replica)
168  ! Read in replicas coordinates
169  skip_vel_section = (input_nr_replica /= neb_nr_replica)
170  IF ((iw > 0) .AND. skip_vel_section) THEN
171  WRITE (iw, '(T2,A)') 'NEB| The number of replica in the input is different from the number', &
172  'NEB| of replica requested for NEB. More Replica will be interpolated.', &
173  'NEB| Therefore the possibly provided velocities will not be read.'
174  END IF
175  ! Further check on velocity section...
176  DO i_rep = 1, input_nr_replica
177  vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", &
178  i_rep_section=i_rep)
179  CALL section_vals_get(vel_section, explicit=explicit)
180  skip_vel_section = skip_vel_section .OR. (.NOT. explicit)
181  END DO
182  ! Setup cartesian coordinates and COLVAR (if requested)
183  coords%xyz(:, :) = 0.0_dp
184  DO i_rep = 1, input_nr_replica
185  coord_section => section_vals_get_subs_vals(replica_section, "COORD", &
186  i_rep_section=i_rep)
187  CALL section_vals_get(coord_section, explicit=explicit)
188  ! Cartesian Coordinates
189  IF (explicit) THEN
190  CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", &
191  n_rep_val=natom)
192  cpassert((natom == SIZE(particle_set)))
193  DO iatom = 1, natom
194  CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", &
195  i_rep_val=iatom, r_vals=rptr)
196  ic = 3*(iatom - 1)
197  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  shell_index = particle_set(iatom)%shell_index
200  IF (shell_index /= 0) THEN
201  is = 3*(natom + shell_index - 1)
202  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  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  i_rep_section=i_rep, c_val=filename)
212  cpassert(trim(filename) /= "")
213  CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.true.)
214  CALL parser_get_next_line(parser, 1)
215  ! Start parser
216  CALL parser_get_object(parser, natom)
217  cpassert((natom == SIZE(particle_set)))
218  CALL parser_get_next_line(parser, 1)
219  DO iatom = 1, natom
220  ! Atom coordinates
221  CALL parser_get_next_line(parser, 1, at_end=my_end)
222  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  " line with title is missing or is empty. Please check the XYZ file and rerun your job!")
227  READ (parser%input_line, *) dummy_char, r(1:3)
228  ic = 3*(iatom - 1)
229  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  shell_index = particle_set(iatom)%shell_index
232  IF (shell_index /= 0) THEN
233  is = 3*(natom + shell_index - 1)
234  coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep)
235  END IF
236  END DO
237  CALL parser_release(parser)
238  END block
239  END IF
240  ! Collective Variables
241  IF (neb_env%use_colvar) THEN
242  CALL section_vals_val_get(replica_section, "COLLECTIVE", &
243  i_rep_section=i_rep, n_rep_val=n_rep)
244  IF (n_rep /= 0) THEN
245  ! Read the values of the collective variables
246  NULLIFY (initial_colvars)
247  CALL section_vals_val_get(replica_section, "COLLECTIVE", &
248  i_rep_section=i_rep, r_vals=initial_colvars)
249  check = (neb_env%nsize_int == SIZE(initial_colvars))
250  cpassert(check)
251  coords%int(:, i_rep) = initial_colvars
252  ELSE
253  ! Compute the values of the collective variables
254  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  CALL dump_replica_coordinates(particle_set, coords, i_rep, i_rep, iw, neb_env%use_colvar)
259  ! Setup Velocities
260  IF (skip_vel_section) THEN
261  CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, &
262  i_rep, iw, globenv, neb_env)
263  ELSE
264  vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", &
265  i_rep_section=i_rep)
266  CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
267  n_rep_val=nval)
268  ! Setup Velocities for collective or cartesian coordinates
269  IF (neb_env%use_colvar) THEN
270  nvar = SIZE(vels%wrk, 1)
271  cpassert(nval == nvar)
272  DO ivar = 1, nvar
273  CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
274  i_rep_val=ivar, r_vals=rptr)
275  vels%wrk(ivar, i_rep) = rptr(1)
276  END DO
277  ELSE
278  natom = SIZE(particle_set)
279  cpassert(nval == natom)
280  DO iatom = 1, natom
281  CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
282  i_rep_val=iatom, r_vals=rptr)
283  ic = 3*(iatom - 1)
284  vels%wrk(ic + 1:ic + 3, i_rep) = rptr(1:3)
285  ! Initially set shell velocities to core velocity
286  shell_index = particle_set(iatom)%shell_index
287  IF (shell_index /= 0) THEN
288  is = 3*(natom + shell_index - 1)
289  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  ALLOCATE (distance(neb_nr_replica - 1))
296  IF (input_nr_replica < neb_nr_replica) THEN
297  ! Interpolate missing replicas
298  nr_replica_to_interpolate = neb_nr_replica - input_nr_replica
299  distance = 0.0_dp
300  IF (iw > 0) THEN
301  WRITE (iw, '(T2,A,I0,A)') 'NEB| Interpolating ', nr_replica_to_interpolate, ' missing Replica.'
302  END IF
303  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  DO j = 1, input_nr_replica - 1
307  CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
308  rotate=neb_env%align_frames)
309  END DO
310  jtarg = maxloc(distance(1:input_nr_replica), 1)
311  IF (iw > 0) THEN
312  WRITE (iw, '(T2,3(A,I0),A)') 'NEB| Interpolating Nr. ', &
313  nr_replica_to_interpolate, ' missing Replica between Replica Nr. ', &
314  jtarg, ' and ', jtarg + 1, '.'
315  END IF
316  input_nr_replica = input_nr_replica + 1
317  nr_replica_to_interpolate = nr_replica_to_interpolate - 1
318  ! Interpolation is a simple bisection in XYZ
319  coords%xyz(:, jtarg + 2:input_nr_replica) = coords%xyz(:, jtarg + 1:input_nr_replica - 1)
320  coords%xyz(:, jtarg + 1) = (coords%xyz(:, jtarg) + coords%xyz(:, jtarg + 2))/2.0_dp
321  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  coords%int(:, jtarg + 2:input_nr_replica) = coords%int(:, jtarg + 1:input_nr_replica - 1)
327  coords%int(:, jtarg + 1) = (coords%int(:, jtarg) + coords%int(:, jtarg + 2))/2.0_dp
328  END IF
329  vels%wrk(:, jtarg + 2:input_nr_replica) = vels%wrk(:, jtarg + 1:input_nr_replica - 1)
330  vels%wrk(:, jtarg + 1) = 0.0_dp
331  CALL dump_replica_coordinates(particle_set, coords, jtarg + 1, &
332  input_nr_replica, iw, neb_env%use_colvar)
333  CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, &
334  jtarg + 1, iw, globenv, neb_env)
335  END DO
336  END IF
337  vels%wrk(:, 1) = 0.0_dp
338  vels%wrk(:, neb_nr_replica) = 0.0_dp
339  ! If we perform a DIIS optimization we don't need velocities
340  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  DO j = 1, input_nr_replica - 1
344  CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
345  rotate=neb_env%align_frames)
346  END DO
347  DEALLOCATE (distance)
348 
349  CALL timestop(handle)
350 
351  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tangent, tmp_a, tmp_b
381  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cvalues, mmatrix, mmatrix_tmp
382 
383  CALL timeset(routinen, handle)
384  n_int = neb_env%nsize_int
385  n_rep_neb = neb_env%number_of_replica
386  n_rep = rep_env%nrep
387  nsize_wrk = coords%size_wrk(1)
388  energies = 0.0_dp
389  ALLOCATE (cvalues(n_int, n_rep))
390  ALLOCATE (mmatrix_tmp(n_int*n_int, n_rep))
391  ALLOCATE (mmatrix(n_int*n_int, n_rep_neb))
392  IF (output_unit > 0) WRITE (output_unit, '(/,T2,A)') "NEB| Computing Energies and Forces"
393  DO irep = 1, n_rep_neb, n_rep
394  DO j = 0, n_rep - 1
395  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  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  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  SELECT CASE (neb_env%pot_type)
408  CASE (pot_neb_full)
409  ! Full potential Energy
410  CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
411  CASE (pot_neb_fe)
412  ! Free Energy Case
413  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  CALL perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, mmatrix_tmp)
417  END SELECT
418 
419  DO j = 0, n_rep - 1
420  IF (irep + j <= n_rep_neb) THEN
421  ! Copy back Forces and Energies
422  forces%wrk(:, irep + j) = rep_env%f(1:nsize_wrk, j + 1)
423  energies(irep + j) = rep_env%f(rep_env%ndim + 1, j + 1)
424  SELECT CASE (neb_env%pot_type)
425  CASE (pot_neb_full)
426  ! Dump Info
427  IF (output_unit > 0) THEN
428  WRITE (output_unit, '(T2,A,I5,A,I5,A)') &
429  "NEB| REPLICA Nr.", irep + j, "- Energy and Forces"
430  WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
431  "NEB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
432  WRITE (output_unit, '(T2,"NEB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
433  DO i = 1, SIZE(particle_set)
434  WRITE (output_unit, '(T2,"NEB|",T12,A,T30,3(2X,F15.9))') &
435  particle_set(i)%atomic_kind%name, &
436  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  coords%xyz(:, irep + j) = rep_env%r(1:rep_env%ndim, j + 1)
443  mmatrix(:, irep + j) = mmatrix_tmp(:, j + 1)
444  IF (output_unit > 0) THEN
445  WRITE (output_unit, '(/,T2,A,I5,A,I5,A)') &
446  "NEB| REPLICA Nr.", irep + j, "- Energy, Collective Variables, Forces"
447  WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
448  "NEB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
449  WRITE (output_unit, &
450  '(T2,"NEB|",T10,"CV Nr.",12X,"Expected COLVAR",5X,"Present COLVAR",10X,"Forces")')
451  DO i = 1, n_int
452  WRITE (output_unit, '(T2,"NEB|",T12,I2,7X,3(5X,F15.9))') &
453  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  DEALLOCATE (cvalues)
461  DEALLOCATE (mmatrix_tmp)
462  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  neb_env%nr_HE_image = maxloc(energies(2:n_rep_neb - 1), 1) + 1
466  ALLOCATE (tangent(nsize_wrk))
467  ! Then modify image forces accordingly to the scheme chosen for the
468  ! calculation.
469  neb_env%spring_energy = 0.0_dp
470  IF (neb_env%optimize_end_points) THEN
471  ALLOCATE (tmp_a(SIZE(forces%wrk, 1)))
472  ALLOCATE (tmp_b(SIZE(forces%wrk, 1)))
473  tmp_a(:) = forces%wrk(:, 1)
474  tmp_b(:) = forces%wrk(:, SIZE(forces%wrk, 2))
475  END IF
476  DO i = 2, neb_env%number_of_replica
477  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  iw=output_unit)
480  END DO
481  IF (neb_env%optimize_end_points) THEN
482  forces%wrk(:, 1) = tmp_a ! Image A
483  forces%wrk(:, SIZE(forces%wrk, 2)) = tmp_b ! Image B
484  DEALLOCATE (tmp_a)
485  DEALLOCATE (tmp_b)
486  ELSE
487  ! Nullify forces on the two end points images
488  forces%wrk(:, 1) = 0.0_dp ! Image A
489  forces%wrk(:, SIZE(forces%wrk, 2)) = 0.0_dp ! Image B
490  END IF
491  DEALLOCATE (tangent)
492  END IF
493  DEALLOCATE (mmatrix)
494  CALL timestop(handle)
495  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  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  CALL timeset(routinen, handle)
524  CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
525  handle=handle2)
526  logger => cp_get_default_logger()
527  CALL force_env_get(f_env%force_env, globenv=globenv, &
528  root_section=root_section)
529  j = rep_env%local_rep_indices(1) - 1
530  n_el = 3*rep_env%nparticle
531  mmatrix = 0.0_dp
532  ! Syncronize position on the replica procs
533  CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
534  cpassert(ierr == 0)
535  !
536  IF (irep + j <= n_rep_neb) THEN
537  logger%iter_info%iteration(2) = irep + j
538  CALL remove_restart_info(root_section)
539  md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
540  CALL section_vals_get(md_section, explicit=explicit)
541  cpassert(explicit)
542  ! Let's syncronize the target of Collective Variables for this run
543  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  CALL qs_mol_dyn(f_env%force_env, globenv=globenv)
547  ! Collect the equilibrated coordinates
548  CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
549  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  cpabort("")
554  rep_env%f(:, j + 1) = 0.0_dp
555  mmatrix = 0.0_dp
556  ELSE
557  rep_env%r(:, j + 1) = 0.0_dp
558  rep_env%f(:, j + 1) = 0.0_dp
559  cvalues(:, j + 1) = 0.0_dp
560  mmatrix(:, j + 1) = 0.0_dp
561  END IF
562  CALL rep_env_sync(rep_env, rep_env%f)
563  CALL rep_env_sync(rep_env, rep_env%r)
564  CALL rep_env_sync(rep_env, cvalues)
565  CALL rep_env_sync(rep_env, mmatrix)
566  CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
567  cpassert(ierr == 0)
568  CALL timestop(handle)
569  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  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  CALL timeset(routinen, handle)
599  CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
600  handle=handle2)
601  logger => cp_get_default_logger()
602  CALL force_env_get(f_env%force_env, globenv=globenv, &
603  root_section=root_section)
604  j = rep_env%local_rep_indices(1) - 1
605  n_el = 3*rep_env%nparticle
606  mmatrix = 0.0_dp
607  ! Syncronize position on the replica procs
608  CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
609  cpassert(ierr == 0)
610  IF (irep + j <= n_rep_neb) THEN
611  logger%iter_info%iteration(2) = irep + j
612  CALL remove_restart_info(root_section)
613  geoopt_section => section_vals_get_subs_vals(root_section, "MOTION%GEO_OPT")
614  CALL section_vals_get(geoopt_section, explicit=explicit)
615  cpassert(explicit)
616  ! Let's syncronize the target of Collective Variables for this run
617  CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env)
618  ! Do a geometry optimization..
619  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  calc_force=.true., skip_external_control=.true.)
624  ! Collect the optimized coordinates
625  CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
626  cpassert(ierr == 0)
627  ! Collect the gradients in cartesian coordinates
628  CALL get_force(rep_env%f_env_id, rep_env%f(1:n_el, j + 1), n_el, ierr)
629  cpassert(ierr == 0)
630  ! Copy the energy
631  CALL get_energy(rep_env%f_env_id, rep_env%f(n_el + 1, j + 1), ierr)
632  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  SIZE(coords%xyz, 1), SIZE(coords%wrk, 1), cvalues(:, j + 1), mmatrix(:, j + 1))
636  ELSE
637  rep_env%r(:, j + 1) = 0.0_dp
638  rep_env%f(:, j + 1) = 0.0_dp
639  cvalues(:, j + 1) = 0.0_dp
640  mmatrix(:, j + 1) = 0.0_dp
641  END IF
642  CALL rep_env_sync(rep_env, rep_env%f)
643  CALL rep_env_sync(rep_env, rep_env%r)
644  CALL rep_env_sync(rep_env, cvalues)
645  CALL rep_env_sync(rep_env, mmatrix)
646  CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
647  cpassert(ierr == 0)
648  CALL timestop(handle)
649  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  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  cpassert(ASSOCIATED(coords))
673  tangent(:) = 0.0_dp
674  ! For the last point we don't need any tangent..
675  IF (i == neb_env%number_of_replica) RETURN
676  ! Several kind of tangents implemented...
677  SELECT CASE (neb_env%id_type)
678  CASE (do_eb)
679  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  rotate=.false.)
683  CALL neb_replica_distance(coords=coords, i0=i + 1, i=i, distance=distance2, iw=iw, &
684  rotate=.false.)
685  tangent(:) = (coords%wrk(:, i) - coords%wrk(:, i - 1))/distance1 + &
686  (coords%wrk(:, i + 1) - coords%wrk(:, i))/distance2
687  CASE (do_it_neb, do_ci_neb, do_d_neb)
688  IF ((energies(i + 1) .GT. energies(i)) .AND. (energies(i) .GT. (energies(i - 1)))) THEN
689  tangent(:) = coords%wrk(:, i + 1) - coords%wrk(:, i)
690  ELSE IF ((energies(i + 1) .LT. energies(i)) .AND. (energies(i) .LT. (energies(i - 1)))) THEN
691  tangent(:) = coords%wrk(:, i) - coords%wrk(:, i - 1)
692  ELSE
693  dvmax = max(abs(energies(i + 1) - energies(i)), abs(energies(i - 1) - energies(i)))
694  dvmin = min(abs(energies(i + 1) - energies(i)), abs(energies(i - 1) - energies(i)))
695  IF (energies(i + 1) .GE. energies(i - 1)) THEN
696  tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*dvmax + (coords%wrk(:, i) - coords%wrk(:, i - 1))*dvmin
697  ELSE
698  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  tangent(:) = 0.0_dp
704  END SELECT
705  distance0 = sqrt(dot_product(tangent(:), tangent(:)))
706  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dtmp1, wrk
735 
736  my_tag = neb_env%id_type
737  IF (PRESENT(tag)) my_tag = tag
738  cpassert(ASSOCIATED(forces))
739  cpassert(ASSOCIATED(coords))
740  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  SELECT CASE (my_tag)
744  CASE (do_b_neb, do_it_neb, do_ci_neb, do_d_neb)
745  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  CALL cite_reference(e2002)
751  RETURN
752  END SELECT
753  ! otherwise proceeed normally..
754  ALLOCATE (wrk(nsize_wrk))
755  ! Spring Energy
756  CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance1, iw=iw, &
757  rotate=.false.)
758  tmp = distance1 - neb_env%avg_distance
759  neb_env%spring_energy = neb_env%spring_energy + 0.5_dp*neb_env%k*tmp**2
760  SELECT CASE (my_tag)
761  CASE (do_eb)
762  CALL cite_reference(elber1987)
763  ! Elastic band - Hamiltonian formulation according the original Karplus/Elber
764  ! formulation
765  ALLOCATE (dtmp1(nsize_wrk))
766  ! derivatives of the spring
767  tmp = distance1 - neb_env%avg_distance
768  dtmp1(:) = 1.0_dp/distance1*(coords%wrk(:, i) - coords%wrk(:, i - 1))
769  wrk(:) = neb_env%k*tmp*dtmp1
770  forces%wrk(:, i) = forces%wrk(:, i) - wrk
771  forces%wrk(:, i - 1) = forces%wrk(:, i - 1) + wrk
772  ! derivatives due to the average length of the spring
773  fac = 1.0_dp/(neb_env%avg_distance*real(neb_env%number_of_replica - 1, kind=dp))
774  wrk(:) = neb_env%k*fac*(coords%wrk(:, i) - coords%wrk(:, i - 1))
775  tmp = 0.0_dp
776  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  rotate=.false.)
779  tmp = tmp + distance1 - neb_env%avg_distance
780  END DO
781  forces%wrk(:, i) = forces%wrk(:, i) + wrk*tmp
782  forces%wrk(:, i - 1) = forces%wrk(:, i - 1) - wrk*tmp
783  DEALLOCATE (dtmp1)
784  CASE (do_b_neb)
785  ! Bisection NEB
786  CALL cite_reference(jonsson1998)
787  wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
788  tmp = neb_env%k*dot_product(wrk, tangent)
789  wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, mmatrix)*tangent
790  forces%wrk(:, i) = wrk + tmp*tangent
791  CASE (do_it_neb)
792  ! Improved tangent NEB
793  CALL cite_reference(jonsson2000_1)
794  CALL neb_replica_distance(coords=coords, i0=i, i=i + 1, distance=distance1, iw=iw, &
795  rotate=.false.)
796  CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance2, iw=iw, &
797  rotate=.false.)
798  tmp = neb_env%k*(distance1 - distance2)
799  wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, mmatrix)*tangent
800  forces%wrk(:, i) = wrk + tmp*tangent
801  CASE (do_ci_neb)
802  ! Climbing Image NEB
803  CALL cite_reference(jonsson2000_2)
804  IF (neb_env%istep <= neb_env%nsteps_it .OR. i /= neb_env%nr_HE_image) THEN
805  CALL get_neb_force(neb_env, tangent, coords, i, forces, do_it_neb, mmatrix, iw)
806  ELSE
807  wrk(:) = forces%wrk(:, i)
808  tmp = -2.0_dp*dot_product_band(neb_env, wrk, tangent, mmatrix)
809  forces%wrk(:, i) = wrk + tmp*tangent
810  END IF
811  CASE (do_d_neb)
812  ! Doubly NEB
813  CALL cite_reference(wales2004)
814  ALLOCATE (dtmp1(nsize_wrk))
815  dtmp1(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, mmatrix)*tangent
816  forces%wrk(:, i) = dtmp1
817  tmp = sqrt(dot_product(dtmp1, dtmp1))
818  dtmp1(:) = dtmp1(:)/tmp
819  ! Project out only the spring component interfering with the
820  ! orthogonal gradient of the band
821  wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
822  tmp = dot_product(wrk, dtmp1)
823  dtmp1(:) = neb_env%k*(wrk(:) - tmp*dtmp1(:))
824  forces%wrk(:, i) = forces%wrk(:, i) + dtmp1(:)
825  DEALLOCATE (dtmp1)
826  END SELECT
827  DEALLOCATE (wrk)
828  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  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  IF (neb_env%use_colvar) THEN
851  nsize_int = neb_env%nsize_int
852  check = ((SIZE(array1) /= SIZE(array2)) .OR. &
853  (SIZE(array1) /= nsize_int) .OR. &
854  (SIZE(array3) /= nsize_int*nsize_int))
855  ! This condition should always be satisfied..
856  cpassert(check)
857  value = dot_product(matmul(reshape(array3, (/nsize_int, nsize_int/)), array1), array2)
858  ELSE
859  value = dot_product(array1, array2)
860  END IF
861  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  SUBROUTINE reorient_images(rotate_frames, particle_set, coords, vels, iw, &
876  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  rot = 0.0_dp
892  rot(1, 1) = 1.0_dp
893  rot(2, 2) = 1.0_dp
894  rot(3, 3) = 1.0_dp
895  DO i = 2, number_of_replica
896  ! The rotation of the replica is enabled exclusively when working in
897  ! cartesian coordinates
898  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  rotate=.true., rot=rot)
901  ! Rotate velocities
902  DO k = 1, SIZE(vels%xyz, 1)/3
903  kind = (k - 1)*3
904  tmp = vels%xyz(kind + 1:kind + 3, i)
905  vels%xyz(kind + 1:kind + 3, i) = matmul(transpose(rot), tmp)
906  END DO
907  END IF
908  IF (PRESENT(distances)) THEN
909  check = SIZE(distances) == (number_of_replica - 1)
910  cpassert(check)
911  xtmp = dot_product(coords%wrk(:, i) - coords%wrk(:, i - 1), &
912  coords%wrk(:, i) - coords%wrk(:, i - 1))
913  distances(i - 1) = sqrt(xtmp)
914  END IF
915  END DO
916  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  SUBROUTINE reparametrize_images(reparametrize_frames, spline_order, smoothing, &
929  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_coords
940 
941  IF (reparametrize_frames) THEN
942  ALLOCATE (tmp_coords(SIZE(coords, 1), SIZE(coords, 2)))
943  tmp_coords(:, :) = coords
944  ! Smoothing
945  DO i = 2, SIZE(coords, 2) - 1
946  coords(:, i) = tmp_coords(:, i)*(1.0_dp - 2.0_dp*smoothing) + &
947  tmp_coords(:, i - 1)*smoothing + tmp_coords(:, i + 1)*smoothing
948  END DO
949  sline = coords - tmp_coords + sline
950  tmp_coords(:, :) = coords
951  ! Reparametrization
952  SELECT CASE (spline_order)
953  CASE (1)
954  ! Compute distances
955  DO i = 2, SIZE(coords, 2)
956  xtmp = dot_product(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
957  distances(i - 1) = sqrt(xtmp)
958  END DO
959  avg_distance = sum(distances)/real(SIZE(coords, 2) - 1, kind=dp)
960  ! Redistribute frames
961  DO i = 2, SIZE(coords, 2) - 1
962  xtmp = 0.0_dp
963  DO j = 1, SIZE(coords, 2) - 1
964  xtmp = xtmp + distances(j)
965  IF (xtmp > avg_distance*real(i - 1, kind=dp)) THEN
966  xtmp = (xtmp - avg_distance*real(i - 1, kind=dp))/distances(j)
967  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  DO i = 2, SIZE(coords, 2)
974  xtmp = dot_product(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
975  distances(i - 1) = sqrt(xtmp)
976  END DO
977  CASE DEFAULT
978  cpwarn("String Method: Spline order greater than 1 not implemented.")
979  END SELECT
980  sline = coords - tmp_coords + sline
981  DEALLOCATE (tmp_coords)
982  END IF
983  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  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  NULLIFY (logger, cc_section)
1007  logger => cp_get_default_logger()
1008  cc_section => section_vals_get_subs_vals(neb_env%neb_section, "CONVERGENCE_CONTROL")
1009  CALL section_vals_val_get(cc_section, "MAX_DR", r_val=max_dr)
1010  CALL section_vals_val_get(cc_section, "MAX_FORCE", r_val=max_force)
1011  CALL section_vals_val_get(cc_section, "RMS_DR", r_val=rms_dr)
1012  CALL section_vals_val_get(cc_section, "RMS_FORCE", r_val=rms_force)
1013  converged = .false.
1014  labels = " NO"
1015  my_max_dr = maxval(abs(dcoords%wrk))
1016  my_max_force = maxval(abs(forces%wrk))
1017  my_rms_dr = sqrt(sum(dcoords%wrk*dcoords%wrk)/real(SIZE(dcoords%wrk, 1)*SIZE(dcoords%wrk, 2), kind=dp))
1018  my_rms_force = sqrt(sum(forces%wrk*forces%wrk)/real(SIZE(forces%wrk, 1)*SIZE(forces%wrk, 2), kind=dp))
1019  IF (my_max_dr < max_dr) labels(1) = "YES"
1020  IF (my_max_force < max_force) labels(2) = "YES"
1021  IF (my_rms_dr < rms_dr) labels(3) = "YES"
1022  IF (my_rms_force < rms_force) labels(4) = "YES"
1023  IF (all(labels == "YES")) converged = .true.
1024 
1025  iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "CONVERGENCE_INFO", &
1026  extension=".nebLog")
1027  IF (iw > 0) THEN
1028  ! Print convergence info
1029  WRITE (iw, fmt='(A,A)') ' **************************************', &
1030  '*****************************************'
1031  WRITE (iw, fmt='(1X,A,2X,F8.5,5X,"[",F8.5,"]",1X,T76,"(",A,")")') &
1032  'RMS DISPLACEMENT =', my_rms_dr, rms_dr, labels(3), &
1033  'MAX DISPLACEMENT =', my_max_dr, max_dr, labels(1), &
1034  'RMS FORCE =', my_rms_force, rms_force, labels(4), &
1035  'MAX FORCE =', my_max_force, max_force, labels(2)
1036  WRITE (iw, fmt='(A,A)') ' **************************************', &
1037  '*****************************************'
1038  END IF
1039  CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, &
1040  "CONVERGENCE_INFO")
1041  END FUNCTION check_convergence
1042 
1043 END MODULE neb_utils
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public e2002
Definition: bibliography.F:43
integer, save, public elber1987
Definition: bibliography.F:43
integer, save, public jonsson2000_1
Definition: bibliography.F:43
integer, save, public jonsson1998
Definition: bibliography.F:43
integer, save, public jonsson2000_2
Definition: bibliography.F:43
integer, save, public wales2004
Definition: bibliography.F:43
evaluations of colvar for internal coordinates schemes
Definition: colvar_utils.F:14
subroutine, public eval_colvar(force_env, coords, cvalues, Bmatrix, MassI, Amatrix)
Computes the values of colvars and the Wilson matrix B and its invers A.
Definition: colvar_utils.F:192
subroutine, public get_clv_force(force_env, forces, coords, nsize_xyz, nsize_int, cvalues, Mmatrix)
Computes the forces in the frame of collective variables, and additional also the local metric tensor...
Definition: colvar_utils.F:512
subroutine, public set_colvars_target(targets, force_env)
Set the value of target for constraints/restraints.
Definition: colvar_utils.F:133
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
interface to use cp2k as library
Definition: f77_interface.F:20
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
subroutine, public get_energy(env_id, e_pot, ierr)
returns the energy of the last configuration calculated
subroutine, public get_pos(env_id, pos, n_el, ierr)
gets the positions of the particles
subroutine, public set_pos(env_id, new_pos, n_el, ierr)
sets the positions of the particles
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
subroutine, public get_force(env_id, frc, n_el, ierr)
gets the forces of the particles
Interface for the force calculations.
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
performs geometry optimization
Definition: geo_opt.F:13
subroutine, public cp_geo_opt(force_env, globenv, eval_opt_geo, rm_restart_info)
Main driver to perform geometry optimization.
Definition: geo_opt.F:58
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public pot_neb_me
integer, parameter, public do_band_cartesian
integer, parameter, public pot_neb_fe
integer, parameter, public do_sm
integer, parameter, public do_b_neb
integer, parameter, public do_d_neb
integer, parameter, public do_eb
integer, parameter, public band_diis_opt
integer, parameter, public do_ci_neb
integer, parameter, public do_it_neb
integer, parameter, public pot_neb_full
checks the input and perform some automatic "magic" on it
subroutine, public remove_restart_info(input_file)
Removes section used to restart a calculation from an input file in memory.
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
integer, parameter, public default_path_length
Definition: kinds.F:58
Perform a molecular dynamics (MD) run using QUICKSTEP.
Definition: md_run.F:14
subroutine, public qs_mol_dyn(force_env, globenv, averages, rm_restart_info, hmc_e_initial, hmc_e_final, mdctrl)
Main driver module for Molecular Dynamics.
Definition: md_run.F:121
Interface to the message passing library MPI.
I/O Module for Nudged Elastic Band Calculation.
Definition: neb_io.F:20
subroutine, public dump_replica_coordinates(particle_set, coords, i_rep, ienum, iw, use_colvar)
dump coordinates of a replica NEB
Definition: neb_io.F:375
subroutine, public handle_band_file_names(rep_env, irep, n_rep, istep)
Handles the correct file names during a band calculation.
Definition: neb_io.F:413
Module with utility to perform MD Nudged Elastic Band Calculation.
Definition: neb_md_utils.F:20
subroutine, public neb_initialize_velocity(vels, neb_section, particle_set, i_rep, iw, globenv, neb_env)
Initialize velocities of replica in an MD optimized algorithm within NEB.
Definition: neb_md_utils.F:66
Typo for Nudged Elastic Band Calculation.
Definition: neb_types.F:20
Module with utility for Nudged Elastic Band Calculation.
Definition: neb_utils.F:20
subroutine, public neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, particle_set, output_unit)
Driver to compute energy and forces within a NEB, Based on the use of the replica_env.
Definition: neb_utils.F:367
logical function, public check_convergence(neb_env, Dcoords, forces)
Checks for convergence criteria during a NEB run.
Definition: neb_utils.F:994
subroutine, public build_replica_coords(neb_section, particle_set, coords, vels, neb_env, iw, globenv, para_env)
Constructs or Read the coordinates for all replica.
Definition: neb_utils.F:140
subroutine, public reorient_images(rotate_frames, particle_set, coords, vels, iw, distances, number_of_replica)
Reorient iteratively all images of the NEB chain in order to have always the smaller RMSD between two...
Definition: neb_utils.F:877
subroutine, public reparametrize_images(reparametrize_frames, spline_order, smoothing, coords, sline, distances)
Reparametrization of the replica for String Method with splines.
Definition: neb_utils.F:930
Define the data structure for the particle information.
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
methods to setup replicas of the same system differing only by atom positions and velocities (as used...
subroutine, public rep_env_calc_e_f(rep_env, calc_f)
evaluates the forces
types used to handle many replica of the same system that differ only in atom positions,...
Definition: replica_types.F:21
subroutine, public rep_env_sync(rep_env, vals)
sends the data from each replica to all the other on replica j/=i data from replica i overwrites val(...
Defines functions to perform rmsd in 3D.
Definition: rmsd.F:12
subroutine, public rmsd3(particle_set, r, r0, output_unit, weights, my_val, rotate, transl, rot, drmsd3)
Computes the RMSD in 3D. Provides also derivatives.
Definition: rmsd.F:53