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