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