(git:ed6f26b)
Loading...
Searching...
No Matches
neb_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
20MODULE neb_utils
21 USE bibliography, ONLY: e2002,&
22 elber1987,&
23 jonsson1998,&
24 jonsson2000_1,&
25 jonsson2000_2,&
26 wales2004,&
27 cite_reference
28 USE colvar_utils, ONLY: eval_colvar,&
29 get_clv_force,&
30 set_colvars_target
31 USE cp_log_handling, ONLY: cp_get_default_logger,&
32 cp_logger_type
33 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
34 cp_print_key_unit_nr
35 USE cp_parser_methods, ONLY: parser_get_next_line,&
36 parser_get_object
37 USE cp_parser_types, ONLY: cp_parser_type,&
38 parser_create,&
39 parser_release
40 USE f77_interface, ONLY: f_env_add_defaults,&
41 f_env_rm_defaults,&
42 f_env_type,&
43 get_energy,&
44 get_force,&
45 get_pos,&
46 set_pos
47 USE force_env_methods, ONLY: force_env_calc_energy_force
48 USE force_env_types, ONLY: force_env_get
49 USE geo_opt, ONLY: cp_geo_opt
50 USE global_types, ONLY: global_environment_type
51 USE input_constants, ONLY: &
52 band_diis_opt, do_b_neb, do_band_cartesian, do_ci_neb, do_d_neb, do_eb, do_it_neb, do_sm, &
53 pot_neb_fe, pot_neb_full, pot_neb_me
54 USE input_cp2k_check, ONLY: remove_restart_info
55 USE input_section_types, ONLY: section_vals_get,&
56 section_vals_get_subs_vals,&
57 section_vals_type,&
58 section_vals_val_get
59 USE kinds, ONLY: default_path_length,&
60 default_string_length,&
61 dp
62 USE md_run, ONLY: qs_mol_dyn
63 USE message_passing, ONLY: mp_para_env_type
64 USE neb_io, ONLY: dump_replica_coordinates,&
65 handle_band_file_names
66 USE neb_md_utils, ONLY: neb_initialize_velocity
67 USE neb_types, ONLY: neb_type,&
68 neb_var_type
69 USE particle_types, ONLY: particle_type
70 USE physcon, ONLY: bohr
71 USE replica_methods, ONLY: rep_env_calc_e_f
72 USE replica_types, ONLY: rep_env_sync,&
73 replica_env_type
74 USE rmsd, ONLY: rmsd3
75#include "../base/base_uses.f90"
76
77 IMPLICIT NONE
78 PRIVATE
79 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_utils'
80 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
81
82 PUBLIC :: build_replica_coords, &
83 neb_calc_energy_forces, &
84 reorient_images, &
85 reparametrize_images, &
86 check_convergence
87
88CONTAINS
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
1043END MODULE neb_utils
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:51
cp_parser_methods::parser_get_object
Definition cp_parser_methods.F:54
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::e2002
integer, save, public e2002
Definition bibliography.F:36
bibliography::elber1987
integer, save, public elber1987
Definition bibliography.F:36
bibliography::jonsson2000_1
integer, save, public jonsson2000_1
Definition bibliography.F:36
bibliography::jonsson1998
integer, save, public jonsson1998
Definition bibliography.F:36
bibliography::jonsson2000_2
integer, save, public jonsson2000_2
Definition bibliography.F:36
bibliography::wales2004
integer, save, public wales2004
Definition bibliography.F:36
colvar_utils
evaluations of colvar for internal coordinates schemes
Definition colvar_utils.F:14
colvar_utils::eval_colvar
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
colvar_utils::get_clv_force
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:513
colvar_utils::set_colvars_target
subroutine, public set_colvars_target(targets, force_env)
Set the value of target for constraints/restraints.
Definition colvar_utils.F:133
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
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)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
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,...
Definition cp_output_handling.F:1063
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
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 ...
Definition cp_parser_methods.F:421
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition cp_parser_types.F:109
cp_parser_types::parser_create
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.
Definition cp_parser_types.F:147
f77_interface
interface to use cp2k as library
Definition f77_interface.F:20
f77_interface::f_env_add_defaults
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 ...
Definition f77_interface.F:452
f77_interface::get_energy
subroutine, public get_energy(env_id, e_pot, ierr)
returns the energy of the last configuration calculated
Definition f77_interface.F:1368
f77_interface::get_pos
subroutine, public get_pos(env_id, pos, n_el, ierr)
gets the positions of the particles
Definition f77_interface.F:1212
f77_interface::set_pos
subroutine, public set_pos(env_id, new_pos, n_el, ierr)
sets the positions of the particles
Definition f77_interface.F:1290
f77_interface::f_env_rm_defaults
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...
Definition f77_interface.F:496
f77_interface::get_force
subroutine, public get_force(env_id, frc, n_el, ierr)
gets the forces of the particles
Definition f77_interface.F:1157
force_env_methods
Interface for the force calculations.
Definition force_env_methods.F:17
force_env_methods::force_env_calc_energy_force
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.
Definition force_env_methods.F:196
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
geo_opt
performs geometry optimization
Definition geo_opt.F:13
geo_opt::cp_geo_opt
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
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::pot_neb_me
integer, parameter, public pot_neb_me
Definition input_constants.F:191
input_constants::do_band_cartesian
integer, parameter, public do_band_cartesian
Definition input_constants.F:67
input_constants::pot_neb_fe
integer, parameter, public pot_neb_fe
Definition input_constants.F:191
input_constants::do_sm
integer, parameter, public do_sm
Definition input_constants.F:184
input_constants::do_b_neb
integer, parameter, public do_b_neb
Definition input_constants.F:184
input_constants::do_d_neb
integer, parameter, public do_d_neb
Definition input_constants.F:184
input_constants::do_eb
integer, parameter, public do_eb
Definition input_constants.F:184
input_constants::band_diis_opt
integer, parameter, public band_diis_opt
Definition input_constants.F:64
input_constants::do_ci_neb
integer, parameter, public do_ci_neb
Definition input_constants.F:184
input_constants::do_it_neb
integer, parameter, public do_it_neb
Definition input_constants.F:184
input_constants::pot_neb_full
integer, parameter, public pot_neb_full
Definition input_constants.F:191
input_cp2k_check
checks the input and perform some automatic "magic" on it
Definition input_cp2k_check.F:14
input_cp2k_check::remove_restart_info
subroutine, public remove_restart_info(input_file)
Removes section used to restart a calculation from an input file in memory.
Definition input_cp2k_check.F:1002
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
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
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
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
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
md_run
Perform a molecular dynamics (MD) run using QUICKSTEP.
Definition md_run.F:14
md_run::qs_mol_dyn
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
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
neb_io
I/O Module for Nudged Elastic Band Calculation.
Definition neb_io.F:20
neb_io::dump_replica_coordinates
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
neb_io::handle_band_file_names
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
neb_md_utils
Module with utility to perform MD Nudged Elastic Band Calculation.
Definition neb_md_utils.F:20
neb_md_utils::neb_initialize_velocity
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
neb_types
Typo for Nudged Elastic Band Calculation.
Definition neb_types.F:20
neb_utils
Module with utility for Nudged Elastic Band Calculation.
Definition neb_utils.F:20
neb_utils::neb_calc_energy_forces
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
neb_utils::check_convergence
logical function, public check_convergence(neb_env, dcoords, forces)
Checks for convergence criteria during a NEB run.
Definition neb_utils.F:994
neb_utils::build_replica_coords
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
neb_utils::reorient_images
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
neb_utils::reparametrize_images
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
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
replica_methods
methods to setup replicas of the same system differing only by atom positions and velocities (as used...
Definition replica_methods.F:16
replica_methods::rep_env_calc_e_f
subroutine, public rep_env_calc_e_f(rep_env, calc_f)
evaluates the forces
Definition replica_methods.F:366
replica_types
types used to handle many replica of the same system that differ only in atom positions,...
Definition replica_types.F:21
replica_types::rep_env_sync
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(...
Definition replica_types.F:229
rmsd
Defines functions to perform rmsd in 3D.
Definition rmsd.F:12
rmsd::rmsd3
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
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_parser_types::cp_parser_type
represent a parser
Definition cp_parser_types.F:77
f77_interface::f_env_type
Definition f77_interface.F:163
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
neb_types::neb_type
Definition neb_types.F:43
neb_types::neb_var_type
Definition neb_types.F:68
particle_types::particle_type
Definition particle_types.F:35
replica_types::replica_env_type
keeps replicated information about the replicas
Definition replica_types.F:96