Line data Source code
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 : ! **************************************************************************************************
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, k, n_rep, natom, &
152 : neb_nr_replica, nr_replica_to_interpolate, nval, nvar, shell_index
153 34 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep_map
154 : LOGICAL :: check, explicit, skip_vel_section
155 34 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: distance
156 : REAL(KIND=dp), DIMENSION(3) :: r
157 34 : REAL(KIND=dp), DIMENSION(:), POINTER :: initial_colvars, rptr
158 : TYPE(section_vals_type), POINTER :: coord_section, replica_section, &
159 : vel_section
160 :
161 34 : CALL timeset(routineN, handle)
162 34 : CPASSERT(ASSOCIATED(coords))
163 34 : CPASSERT(ASSOCIATED(vels))
164 34 : neb_nr_replica = neb_env%number_of_replica
165 34 : replica_section => section_vals_get_subs_vals(neb_section, "REPLICA")
166 34 : CALL section_vals_get(replica_section, n_repetition=input_nr_replica)
167 : ! Calculation is aborted if input replicas are more then the requested ones for the BAND..
168 34 : CPASSERT(input_nr_replica <= neb_nr_replica)
169 : ! Read in replicas coordinates
170 34 : skip_vel_section = (input_nr_replica /= neb_nr_replica)
171 34 : IF ((iw > 0) .AND. skip_vel_section) THEN
172 2 : WRITE (iw, '(T2,A)') 'NEB| The number of replica in the input is different from the number', &
173 2 : 'NEB| of replica requested for NEB. More Replica will be interpolated.', &
174 4 : 'NEB| Therefore the possibly provided velocities will not be read.'
175 : END IF
176 : ! Further check on velocity section...
177 194 : DO i_rep = 1, input_nr_replica
178 : vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", &
179 160 : i_rep_section=i_rep)
180 160 : CALL section_vals_get(vel_section, explicit=explicit)
181 222 : skip_vel_section = skip_vel_section .OR. (.NOT. explicit)
182 : END DO
183 : ! Setup cartesian coordinates and COLVAR (if requested)
184 41772 : coords%xyz(:, :) = 0.0_dp
185 : ! Mapping between input replica and actual replica
186 102 : ALLOCATE (rep_map(neb_nr_replica))
187 34 : rep_map(:) = 0
188 194 : DO i_rep = 1, input_nr_replica
189 : coord_section => section_vals_get_subs_vals(replica_section, "COORD", &
190 160 : i_rep_section=i_rep)
191 160 : CALL section_vals_get(coord_section, explicit=explicit)
192 160 : rep_map(i_rep) = i_rep
193 : ! Cartesian Coordinates
194 160 : IF (explicit) THEN
195 : CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", &
196 32 : n_rep_val=natom)
197 32 : CPASSERT((natom == SIZE(particle_set)))
198 1512 : DO iatom = 1, natom
199 : CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", &
200 1480 : i_rep_val=iatom, r_vals=rptr)
201 1480 : ic = 3*(iatom - 1)
202 10360 : coords%xyz(ic + 1:ic + 3, i_rep) = rptr(1:3)*bohr
203 : ! Initially core and shell positions are set to the atomic positions
204 1480 : shell_index = particle_set(iatom)%shell_index
205 1512 : IF (shell_index /= 0) THEN
206 1140 : is = 3*(natom + shell_index - 1)
207 7980 : coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep)
208 : END IF
209 : END DO
210 : ELSE
211 128 : BLOCK
212 : LOGICAL :: my_end
213 : CHARACTER(LEN=default_string_length) :: dummy_char
214 : TYPE(cp_parser_type) :: parser
215 : CALL section_vals_val_get(replica_section, "COORD_FILE_NAME", &
216 128 : i_rep_section=i_rep, c_val=filename)
217 128 : CPASSERT(TRIM(filename) /= "")
218 128 : CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.TRUE.)
219 128 : CALL parser_get_next_line(parser, 1)
220 : ! Start parser
221 128 : CALL parser_get_object(parser, natom)
222 128 : CPASSERT((natom == SIZE(particle_set)))
223 128 : CALL parser_get_next_line(parser, 1)
224 9774 : DO iatom = 1, natom
225 : ! Atom coordinates
226 9646 : CALL parser_get_next_line(parser, 1, at_end=my_end)
227 9646 : IF (my_end) &
228 : CALL cp_abort(__LOCATION__, &
229 : "Number of lines in XYZ format not equal to the number of atoms."// &
230 : " Error in XYZ format for REPLICA coordinates. Very probably the"// &
231 0 : " line with title is missing or is empty. Please check the XYZ file and rerun your job!")
232 9646 : READ (parser%input_line, *) dummy_char, r(1:3)
233 9646 : ic = 3*(iatom - 1)
234 38584 : 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 9646 : shell_index = particle_set(iatom)%shell_index
237 9774 : IF (shell_index /= 0) THEN
238 0 : is = 3*(natom + shell_index - 1)
239 0 : coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep)
240 : END IF
241 : END DO
242 512 : CALL parser_release(parser)
243 : END BLOCK
244 : END IF
245 : ! Collective Variables
246 160 : IF (neb_env%use_colvar) THEN
247 : CALL section_vals_val_get(replica_section, "COLLECTIVE", &
248 18 : i_rep_section=i_rep, n_rep_val=n_rep)
249 18 : IF (n_rep /= 0) THEN
250 : ! Read the values of the collective variables
251 10 : NULLIFY (initial_colvars)
252 : CALL section_vals_val_get(replica_section, "COLLECTIVE", &
253 10 : i_rep_section=i_rep, r_vals=initial_colvars)
254 10 : check = (neb_env%nsize_int == SIZE(initial_colvars))
255 10 : CPASSERT(check)
256 30 : coords%int(:, i_rep) = initial_colvars
257 : ELSE
258 : ! Compute the values of the collective variables
259 8 : 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 160 : CALL dump_replica_coordinates(particle_set, coords, i_rep, i_rep, iw, neb_env%use_colvar)
264 : ! Setup Velocities
265 354 : IF (skip_vel_section) THEN
266 : CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, &
267 132 : i_rep, iw, globenv, neb_env)
268 : ELSE
269 : vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", &
270 28 : i_rep_section=i_rep)
271 : CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
272 28 : n_rep_val=nval)
273 : ! Setup Velocities for collective or cartesian coordinates
274 28 : IF (neb_env%use_colvar) THEN
275 10 : nvar = SIZE(vels%wrk, 1)
276 10 : CPASSERT(nval == nvar)
277 20 : DO ivar = 1, nvar
278 : CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
279 10 : i_rep_val=ivar, r_vals=rptr)
280 20 : vels%wrk(ivar, i_rep) = rptr(1)
281 : END DO
282 : ELSE
283 18 : natom = SIZE(particle_set)
284 18 : CPASSERT(nval == natom)
285 948 : DO iatom = 1, natom
286 : CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", &
287 930 : i_rep_val=iatom, r_vals=rptr)
288 930 : ic = 3*(iatom - 1)
289 6510 : vels%wrk(ic + 1:ic + 3, i_rep) = rptr(1:3)
290 : ! Initially set shell velocities to core velocity
291 930 : shell_index = particle_set(iatom)%shell_index
292 948 : IF (shell_index /= 0) THEN
293 760 : is = 3*(natom + shell_index - 1)
294 5320 : 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 102 : ALLOCATE (distance(neb_nr_replica - 1))
301 34 : IF (iw > 0) THEN
302 4 : WRITE (iw, '(T2,A)') 'NEB| Mapping between input and requested replica so far'
303 4 : WRITE (iw, '(T2,A)') 'NEB| 1, 2, ... = input replica in sections order,'
304 4 : WRITE (iw, '(T2,A)') 'NEB| -1, -2, ... = added replica in insertion order (if any),'
305 4 : WRITE (iw, '(T2,A)') 'NEB| 0 = yet to be completed replica (if any)'
306 8 : DO j = 1, neb_nr_replica, 8
307 8 : 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 34 : IF (input_nr_replica < neb_nr_replica) THEN
311 : ! Interpolate missing replicas
312 12 : nr_replica_to_interpolate = neb_nr_replica - input_nr_replica
313 12 : k = 0
314 12 : distance = 0.0_dp
315 12 : IF (iw > 0) THEN
316 2 : WRITE (iw, '(T2,A,I0,A)') 'NEB| Interpolating ', nr_replica_to_interpolate, ' missing replica '// &
317 4 : 'by stepwise bisection of distance.'
318 : END IF
319 64 : 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 358 : DO j = 1, input_nr_replica - 1
323 : CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
324 358 : rotate=neb_env%align_frames)
325 : END DO
326 410 : jtarg = MAXLOC(distance(1:input_nr_replica), 1)
327 52 : IF (iw > 0) THEN
328 3 : WRITE (iw, '(/,T2,3(A,I0),A)') 'NEB| Interpolating Nr. ', &
329 3 : nr_replica_to_interpolate, ' missing Replica; next between Replica Nr. ', &
330 6 : jtarg, ' and ', jtarg + 1, '.'
331 : END IF
332 52 : input_nr_replica = input_nr_replica + 1
333 52 : nr_replica_to_interpolate = nr_replica_to_interpolate - 1
334 52 : k = k + 1
335 234 : rep_map(jtarg + 2:input_nr_replica) = rep_map(jtarg + 1:input_nr_replica - 1)
336 52 : IF (jtarg + 1 <= neb_nr_replica .AND. jtarg + 1 >= 1) rep_map(jtarg + 1) = -k
337 : ! Interpolation is a simple bisection in XYZ
338 12630 : coords%xyz(:, jtarg + 2:input_nr_replica) = coords%xyz(:, jtarg + 1:input_nr_replica - 1)
339 4780 : coords%xyz(:, jtarg + 1) = (coords%xyz(:, jtarg) + coords%xyz(:, jtarg + 2))/2.0_dp
340 52 : 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 6 : coords%int(:, jtarg + 2:input_nr_replica) = coords%int(:, jtarg + 1:input_nr_replica - 1)
346 4 : coords%int(:, jtarg + 1) = (coords%int(:, jtarg) + coords%int(:, jtarg + 2))/2.0_dp
347 : END IF
348 12530 : vels%wrk(:, jtarg + 2:input_nr_replica) = vels%wrk(:, jtarg + 1:input_nr_replica - 1)
349 4680 : vels%wrk(:, jtarg + 1) = 0.0_dp
350 : CALL dump_replica_coordinates(particle_set, coords, jtarg + 1, &
351 52 : input_nr_replica, iw, neb_env%use_colvar)
352 : CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, &
353 52 : jtarg + 1, iw, globenv, neb_env)
354 64 : IF (iw > 0) THEN
355 3 : WRITE (iw, '(T2,A)') 'NEB| Mapping between input and requested replica so far'
356 6 : DO j = 1, neb_nr_replica, 8
357 6 : 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 8126 : vels%wrk(:, 1) = 0.0_dp
363 8126 : vels%wrk(:, neb_nr_replica) = 0.0_dp
364 : ! If we perform a DIIS optimization we don't need velocities
365 37092 : 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 212 : DO j = 1, input_nr_replica - 1
369 : CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
370 212 : rotate=neb_env%align_frames)
371 : END DO
372 : ! If there are still large distances, hint that more replica may be requested
373 390 : IF (MAXVAL(distance)/MINVAL(distance) > 1.5_dp) THEN
374 32 : IF (iw > 0) THEN
375 4 : WRITE (iw, '(/,T2,A)') 'NEB| After interpolating replica, the maximum distance between adjacent replica'
376 4 : WRITE (iw, '(T2,A)') 'NEB| is still larger than 1.5 times the minimum distance. In order to make the'
377 4 : WRITE (iw, '(T2,A)') 'NEB| distribution of replica more uniform, consider providing more structures'
378 4 : WRITE (iw, '(T2,A)') 'NEB| with &BAND/&REPLICA sections or increasing &BAND/NUMBER_OF_REPLICA value.'
379 : END IF
380 : END IF
381 34 : DEALLOCATE (distance)
382 34 : DEALLOCATE (rep_map)
383 :
384 34 : CALL timestop(handle)
385 :
386 68 : 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 908 : 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 908 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tangent, tmp_a, tmp_b
416 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cvalues, Mmatrix, Mmatrix_tmp
417 :
418 908 : CALL timeset(routineN, handle)
419 908 : n_int = neb_env%nsize_int
420 908 : n_rep_neb = neb_env%number_of_replica
421 908 : n_rep = rep_env%nrep
422 908 : nsize_wrk = coords%size_wrk(1)
423 5948 : energies = 0.0_dp
424 2748 : ALLOCATE (cvalues(n_int, n_rep))
425 2748 : ALLOCATE (Mmatrix_tmp(n_int*n_int, n_rep))
426 2748 : ALLOCATE (Mmatrix(n_int*n_int, n_rep_neb))
427 908 : IF (output_unit > 0) WRITE (output_unit, '(/,T2,A)') "NEB| Computing Energies and Forces"
428 3700 : DO irep = 1, n_rep_neb, n_rep
429 8346 : DO j = 0, n_rep - 1
430 8346 : 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 2178252 : 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 2792 : 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 5512 : SELECT CASE (neb_env%pot_type)
443 : CASE (pot_neb_full)
444 : ! Full potential Energy
445 2720 : CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
446 : CASE (pot_neb_fe)
447 : ! Free Energy Case
448 0 : 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 2792 : CALL perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix_tmp)
452 : END SELECT
453 :
454 9254 : DO j = 0, n_rep - 1
455 8346 : IF (irep + j <= n_rep_neb) THEN
456 : ! Copy back Forces and Energies
457 2166252 : forces%wrk(:, irep + j) = rep_env%f(1:nsize_wrk, j + 1)
458 5040 : energies(irep + j) = rep_env%f(rep_env%ndim + 1, j + 1)
459 9960 : SELECT CASE (neb_env%pot_type)
460 : CASE (pot_neb_full)
461 : ! Dump Info
462 4920 : IF (output_unit > 0) THEN
463 : WRITE (output_unit, '(T2,A,I5,A,I5,A)') &
464 0 : "NEB| REPLICA Nr.", irep + j, "- Energy and Forces"
465 : WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
466 0 : "NEB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
467 0 : WRITE (output_unit, '(T2,"NEB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
468 0 : DO i = 1, SIZE(particle_set)
469 : WRITE (output_unit, '(T2,"NEB|",T12,A,T30,3(2X,F15.9))') &
470 0 : particle_set(i)%atomic_kind%name, &
471 0 : 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 12480 : coords%xyz(:, irep + j) = rep_env%r(1:rep_env%ndim, j + 1)
478 240 : Mmatrix(:, irep + j) = Mmatrix_tmp(:, j + 1)
479 5160 : IF (output_unit > 0) THEN
480 : WRITE (output_unit, '(/,T2,A,I5,A,I5,A)') &
481 60 : "NEB| REPLICA Nr.", irep + j, "- Energy, Collective Variables, Forces"
482 : WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
483 60 : "NEB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
484 : WRITE (output_unit, &
485 60 : '(T2,"NEB|",T10,"CV Nr.",12X,"Expected COLVAR",5X,"Present COLVAR",10X,"Forces")')
486 120 : DO i = 1, n_int
487 : WRITE (output_unit, '(T2,"NEB|",T12,I2,7X,3(5X,F15.9))') &
488 120 : 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 908 : DEALLOCATE (cvalues)
496 908 : DEALLOCATE (Mmatrix_tmp)
497 908 : 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 4132 : neb_env%nr_HE_image = MAXLOC(energies(2:n_rep_neb - 1), 1) + 1
501 2724 : ALLOCATE (tangent(nsize_wrk))
502 : ! Then modify image forces accordingly to the scheme chosen for the
503 : ! calculation.
504 908 : neb_env%spring_energy = 0.0_dp
505 908 : IF (neb_env%optimize_end_points) THEN
506 1290 : ALLOCATE (tmp_a(SIZE(forces%wrk, 1)))
507 860 : ALLOCATE (tmp_b(SIZE(forces%wrk, 1)))
508 212314 : tmp_a(:) = forces%wrk(:, 1)
509 212314 : tmp_b(:) = forces%wrk(:, SIZE(forces%wrk, 2))
510 : END IF
511 5040 : DO i = 2, neb_env%number_of_replica
512 4132 : 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 5040 : iw=output_unit)
515 : END DO
516 908 : IF (neb_env%optimize_end_points) THEN
517 212314 : forces%wrk(:, 1) = tmp_a ! Image A
518 212314 : forces%wrk(:, SIZE(forces%wrk, 2)) = tmp_b ! Image B
519 430 : DEALLOCATE (tmp_a)
520 430 : DEALLOCATE (tmp_b)
521 : ELSE
522 : ! Nullify forces on the two end points images
523 37102 : forces%wrk(:, 1) = 0.0_dp ! Image A
524 37102 : forces%wrk(:, SIZE(forces%wrk, 2)) = 0.0_dp ! Image B
525 : END IF
526 908 : DEALLOCATE (tangent)
527 : END IF
528 908 : DEALLOCATE (Mmatrix)
529 908 : CALL timestop(handle)
530 1816 : 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 0 : 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 0 : CALL timeset(routineN, handle)
559 : CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
560 0 : handle=handle2)
561 0 : logger => cp_get_default_logger()
562 : CALL force_env_get(f_env%force_env, globenv=globenv, &
563 0 : root_section=root_section)
564 0 : j = rep_env%local_rep_indices(1) - 1
565 0 : n_el = 3*rep_env%nparticle
566 0 : Mmatrix = 0.0_dp
567 : ! Syncronize position on the replica procs
568 0 : CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
569 0 : CPASSERT(ierr == 0)
570 : !
571 0 : IF (irep + j <= n_rep_neb) THEN
572 0 : logger%iter_info%iteration(2) = irep + j
573 0 : CALL remove_restart_info(root_section)
574 0 : md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
575 0 : CALL section_vals_get(md_section, explicit=explicit)
576 0 : CPASSERT(explicit)
577 : ! Let's syncronize the target of Collective Variables for this run
578 0 : 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 0 : CALL qs_mol_dyn(f_env%force_env, globenv=globenv)
582 : ! Collect the equilibrated coordinates
583 0 : CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
584 0 : 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 0 : CPABORT("")
589 0 : rep_env%f(:, j + 1) = 0.0_dp
590 0 : Mmatrix = 0.0_dp
591 : ELSE
592 0 : rep_env%r(:, j + 1) = 0.0_dp
593 0 : rep_env%f(:, j + 1) = 0.0_dp
594 0 : cvalues(:, j + 1) = 0.0_dp
595 0 : Mmatrix(:, j + 1) = 0.0_dp
596 : END IF
597 0 : CALL rep_env_sync(rep_env, rep_env%f)
598 0 : CALL rep_env_sync(rep_env, rep_env%r)
599 0 : CALL rep_env_sync(rep_env, cvalues)
600 0 : CALL rep_env_sync(rep_env, Mmatrix)
601 0 : CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
602 0 : CPASSERT(ierr == 0)
603 0 : CALL timestop(handle)
604 0 : 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 72 : 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 72 : CALL timeset(routineN, handle)
634 : CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
635 72 : handle=handle2)
636 72 : logger => cp_get_default_logger()
637 : CALL force_env_get(f_env%force_env, globenv=globenv, &
638 72 : root_section=root_section)
639 72 : j = rep_env%local_rep_indices(1) - 1
640 72 : n_el = 3*rep_env%nparticle
641 360 : Mmatrix = 0.0_dp
642 : ! Syncronize position on the replica procs
643 72 : CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
644 72 : CPASSERT(ierr == 0)
645 72 : IF (irep + j <= n_rep_neb) THEN
646 60 : logger%iter_info%iteration(2) = irep + j
647 60 : CALL remove_restart_info(root_section)
648 60 : geoopt_section => section_vals_get_subs_vals(root_section, "MOTION%GEO_OPT")
649 60 : CALL section_vals_get(geoopt_section, explicit=explicit)
650 60 : CPASSERT(explicit)
651 : ! Let's syncronize the target of Collective Variables for this run
652 60 : CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env)
653 : ! Do a geometry optimization..
654 60 : 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 60 : calc_force=.TRUE., skip_external_control=.TRUE.)
659 : ! Collect the optimized coordinates
660 60 : CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
661 60 : CPASSERT(ierr == 0)
662 : ! Collect the gradients in cartesian coordinates
663 60 : CALL get_force(rep_env%f_env_id, rep_env%f(1:n_el, j + 1), n_el, ierr)
664 60 : CPASSERT(ierr == 0)
665 : ! Copy the energy
666 60 : CALL get_energy(rep_env%f_env_id, rep_env%f(n_el + 1, j + 1), ierr)
667 60 : 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 60 : SIZE(coords%xyz, 1), SIZE(coords%wrk, 1), cvalues(:, j + 1), Mmatrix(:, j + 1))
671 : ELSE
672 624 : rep_env%r(:, j + 1) = 0.0_dp
673 636 : rep_env%f(:, j + 1) = 0.0_dp
674 24 : cvalues(:, j + 1) = 0.0_dp
675 24 : Mmatrix(:, j + 1) = 0.0_dp
676 : END IF
677 72 : CALL rep_env_sync(rep_env, rep_env%f)
678 72 : CALL rep_env_sync(rep_env, rep_env%r)
679 72 : CALL rep_env_sync(rep_env, cvalues)
680 72 : CALL rep_env_sync(rep_env, Mmatrix)
681 72 : CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
682 72 : CPASSERT(ierr == 0)
683 72 : CALL timestop(handle)
684 72 : 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 4132 : 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 4132 : CPASSERT(ASSOCIATED(coords))
708 833710 : tangent(:) = 0.0_dp
709 : ! For the last point we don't need any tangent..
710 4132 : IF (i == neb_env%number_of_replica) RETURN
711 : ! Several kind of tangents implemented...
712 3224 : SELECT CASE (neb_env%id_type)
713 : CASE (do_eb)
714 38792 : 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 126 : rotate=.FALSE.)
718 : CALL neb_replica_distance(coords=coords, i0=i + 1, i=i, distance=distance2, iw=iw, &
719 126 : rotate=.FALSE.)
720 : tangent(:) = (coords%wrk(:, i) - coords%wrk(:, i - 1))/distance1 + &
721 6552 : (coords%wrk(:, i + 1) - coords%wrk(:, i))/distance2
722 : CASE (do_it_neb, do_ci_neb, do_d_neb)
723 1974 : IF ((energies(i + 1) > energies(i)) .AND. (energies(i) > (energies(i - 1)))) THEN
724 226914 : tangent(:) = coords%wrk(:, i + 1) - coords%wrk(:, i)
725 1282 : ELSE IF ((energies(i + 1) < energies(i)) .AND. (energies(i) < (energies(i - 1)))) THEN
726 35162 : tangent(:) = coords%wrk(:, i) - coords%wrk(:, i - 1)
727 : ELSE
728 1188 : DVmax = MAX(ABS(energies(i + 1) - energies(i)), ABS(energies(i - 1) - energies(i)))
729 1188 : DVmin = MIN(ABS(energies(i + 1) - energies(i)), ABS(energies(i - 1) - energies(i)))
730 1188 : IF (energies(i + 1) >= energies(i - 1)) THEN
731 26028 : tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*DVmax + (coords%wrk(:, i) - coords%wrk(:, i - 1))*DVmin
732 : ELSE
733 231190 : 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 22502 : tangent(:) = 0.0_dp
739 : END SELECT
740 584294 : distance0 = SQRT(DOT_PRODUCT(tangent(:), tangent(:)))
741 526970 : 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 4630 : 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 4630 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dtmp1, wrk
770 :
771 4630 : my_tag = neb_env%id_type
772 4630 : IF (PRESENT(tag)) my_tag = tag
773 4630 : CPASSERT(ASSOCIATED(forces))
774 4630 : CPASSERT(ASSOCIATED(coords))
775 4630 : 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 3336 : SELECT CASE (my_tag)
779 : CASE (do_b_neb, do_it_neb, do_ci_neb, do_d_neb)
780 3336 : 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 420 : CALL cite_reference(E2002)
786 4630 : RETURN
787 : END SELECT
788 : ! otherwise proceeed normally..
789 10416 : ALLOCATE (wrk(nsize_wrk))
790 : ! Spring Energy
791 : CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance1, iw=iw, &
792 3472 : rotate=.FALSE.)
793 3472 : tmp = distance1 - neb_env%avg_distance
794 3472 : neb_env%spring_energy = neb_env%spring_energy + 0.5_dp*neb_env%k*tmp**2
795 874 : SELECT CASE (my_tag)
796 : CASE (do_eb)
797 874 : CALL cite_reference(Elber1987)
798 : ! Elastic band - Hamiltonian formulation according the original Karplus/Elber
799 : ! formulation
800 1748 : ALLOCATE (dtmp1(nsize_wrk))
801 : ! derivatives of the spring
802 874 : tmp = distance1 - neb_env%avg_distance
803 45448 : dtmp1(:) = 1.0_dp/distance1*(coords%wrk(:, i) - coords%wrk(:, i - 1))
804 45448 : wrk(:) = neb_env%k*tmp*dtmp1
805 45448 : forces%wrk(:, i) = forces%wrk(:, i) - wrk
806 45448 : forces%wrk(:, i - 1) = forces%wrk(:, i - 1) + wrk
807 : ! derivatives due to the average length of the spring
808 874 : fac = 1.0_dp/(neb_env%avg_distance*REAL(neb_env%number_of_replica - 1, KIND=dp))
809 45448 : wrk(:) = neb_env%k*fac*(coords%wrk(:, i) - coords%wrk(:, i - 1))
810 874 : tmp = 0.0_dp
811 7880 : 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 7006 : rotate=.FALSE.)
814 7880 : tmp = tmp + distance1 - neb_env%avg_distance
815 : END DO
816 45448 : forces%wrk(:, i) = forces%wrk(:, i) + wrk*tmp
817 45448 : forces%wrk(:, i - 1) = forces%wrk(:, i - 1) - wrk*tmp
818 874 : DEALLOCATE (dtmp1)
819 : CASE (do_b_neb)
820 : ! Bisection NEB
821 126 : CALL cite_reference(Jonsson1998)
822 6552 : wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
823 6552 : tmp = neb_env%k*DOT_PRODUCT(wrk, tangent)
824 6552 : wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent
825 6552 : forces%wrk(:, i) = wrk + tmp*tangent
826 : CASE (do_it_neb)
827 : ! Improved tangent NEB
828 1236 : CALL cite_reference(Jonsson2000_1)
829 : CALL neb_replica_distance(coords=coords, i0=i, i=i + 1, distance=distance1, iw=iw, &
830 1236 : rotate=.FALSE.)
831 : CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance2, iw=iw, &
832 1236 : rotate=.FALSE.)
833 1236 : tmp = neb_env%k*(distance1 - distance2)
834 309648 : wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent
835 309648 : forces%wrk(:, i) = wrk + tmp*tangent
836 : CASE (do_ci_neb)
837 : ! Climbing Image NEB
838 858 : CALL cite_reference(Jonsson2000_2)
839 858 : IF (neb_env%istep <= neb_env%nsteps_it .OR. i /= neb_env%nr_HE_image) THEN
840 498 : CALL get_neb_force(neb_env, tangent, coords, i, forces, do_it_neb, Mmatrix, iw)
841 : ELSE
842 189990 : wrk(:) = forces%wrk(:, i)
843 360 : tmp = -2.0_dp*dot_product_band(neb_env, wrk, tangent, Mmatrix)
844 189990 : forces%wrk(:, i) = wrk + tmp*tangent
845 : END IF
846 : CASE (do_d_neb)
847 : ! Doubly NEB
848 378 : CALL cite_reference(Wales2004)
849 756 : ALLOCATE (dtmp1(nsize_wrk))
850 19656 : dtmp1(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent
851 19656 : forces%wrk(:, i) = dtmp1
852 19656 : tmp = SQRT(DOT_PRODUCT(dtmp1, dtmp1))
853 19656 : dtmp1(:) = dtmp1(:)/tmp
854 : ! Project out only the spring component interfering with the
855 : ! orthogonal gradient of the band
856 19656 : wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
857 19656 : tmp = DOT_PRODUCT(wrk, dtmp1)
858 19656 : dtmp1(:) = neb_env%k*(wrk(:) - tmp*dtmp1(:))
859 19656 : forces%wrk(:, i) = forces%wrk(:, i) + dtmp1(:)
860 3850 : DEALLOCATE (dtmp1)
861 : END SELECT
862 3472 : 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 2100 : 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 2100 : IF (neb_env%use_colvar) THEN
886 72 : nsize_int = neb_env%nsize_int
887 : check = ((SIZE(array1) /= SIZE(array2)) .OR. &
888 : (SIZE(array1) /= nsize_int) .OR. &
889 216 : (SIZE(array3) /= nsize_int*nsize_int))
890 : ! This condition should always be satisfied..
891 0 : CPASSERT(check)
892 576 : value = DOT_PRODUCT(MATMUL(RESHAPE(array3, [nsize_int, nsize_int]), array1), array2)
893 : ELSE
894 525702 : value = DOT_PRODUCT(array1, array2)
895 : END IF
896 2100 : 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 908 : SUBROUTINE reorient_images(rotate_frames, particle_set, coords, vels, iw, &
911 908 : 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 908 : rot = 0.0_dp
927 908 : rot(1, 1) = 1.0_dp
928 908 : rot(2, 2) = 1.0_dp
929 908 : rot(3, 3) = 1.0_dp
930 5040 : DO i = 2, number_of_replica
931 : ! The rotation of the replica is enabled exclusively when working in
932 : ! cartesian coordinates
933 4132 : 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 540 : rotate=.TRUE., rot=rot)
936 : ! Rotate velocities
937 11796 : DO k = 1, SIZE(vels%xyz, 1)/3
938 11256 : kind = (k - 1)*3
939 45024 : tmp = vels%xyz(kind + 1:kind + 3, i)
940 45564 : vels%xyz(kind + 1:kind + 3, i) = MATMUL(TRANSPOSE(rot), tmp)
941 : END DO
942 : END IF
943 5040 : IF (PRESENT(distances)) THEN
944 4132 : check = SIZE(distances) == (number_of_replica - 1)
945 4132 : CPASSERT(check)
946 : xtmp = DOT_PRODUCT(coords%wrk(:, i) - coords%wrk(:, i - 1), &
947 833710 : coords%wrk(:, i) - coords%wrk(:, i - 1))
948 4132 : distances(i - 1) = SQRT(xtmp)
949 : END IF
950 : END DO
951 908 : 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 402 : SUBROUTINE reparametrize_images(reparametrize_frames, spline_order, smoothing, &
964 402 : 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 402 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_coords
975 :
976 402 : IF (reparametrize_frames) THEN
977 168 : ALLOCATE (tmp_coords(SIZE(coords, 1), SIZE(coords, 2)))
978 24066 : tmp_coords(:, :) = coords
979 : ! Smoothing
980 420 : DO i = 2, SIZE(coords, 2) - 1
981 : coords(:, i) = tmp_coords(:, i)*(1.0_dp - 2.0_dp*smoothing) + &
982 19698 : tmp_coords(:, i - 1)*smoothing + tmp_coords(:, i + 1)*smoothing
983 : END DO
984 48090 : sline = coords - tmp_coords + sline
985 24066 : tmp_coords(:, :) = coords
986 : ! Reparametrization
987 84 : SELECT CASE (spline_order)
988 : CASE (1)
989 : ! Compute distances
990 462 : DO i = 2, SIZE(coords, 2)
991 21840 : xtmp = DOT_PRODUCT(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
992 462 : distances(i - 1) = SQRT(xtmp)
993 : END DO
994 462 : avg_distance = SUM(distances)/REAL(SIZE(coords, 2) - 1, KIND=dp)
995 : ! Redistribute frames
996 420 : DO i = 2, SIZE(coords, 2) - 1
997 378 : xtmp = 0.0_dp
998 2022 : DO j = 1, SIZE(coords, 2) - 1
999 1980 : xtmp = xtmp + distances(j)
1000 1980 : IF (xtmp > avg_distance*REAL(i - 1, KIND=dp)) THEN
1001 378 : xtmp = (xtmp - avg_distance*REAL(i - 1, KIND=dp))/distances(j)
1002 19656 : 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 462 : DO i = 2, SIZE(coords, 2)
1009 21840 : xtmp = DOT_PRODUCT(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
1010 462 : distances(i - 1) = SQRT(xtmp)
1011 : END DO
1012 : CASE DEFAULT
1013 42 : CPWARN("String Method: Spline order greater than 1 not implemented.")
1014 : END SELECT
1015 48090 : sline = coords - tmp_coords + sline
1016 42 : DEALLOCATE (tmp_coords)
1017 : END IF
1018 402 : 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 2720 : 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 544 : NULLIFY (cc_section)
1043 544 : cc_section => section_vals_get_subs_vals(neb_env%neb_section, "CONVERGENCE_CONTROL")
1044 544 : CALL section_vals_val_get(cc_section, "MAX_DR", r_val=max_dr)
1045 544 : CALL section_vals_val_get(cc_section, "MAX_FORCE", r_val=max_force)
1046 544 : CALL section_vals_val_get(cc_section, "RMS_DR", r_val=rms_dr)
1047 544 : CALL section_vals_val_get(cc_section, "RMS_FORCE", r_val=rms_force)
1048 544 : converged = .FALSE.
1049 2720 : labels = " NO"
1050 562306 : my_max_dr = MAXVAL(ABS(Dcoords%wrk))
1051 562306 : my_max_force = MAXVAL(ABS(forces%wrk))
1052 562306 : my_rms_dr = SQRT(SUM(Dcoords%wrk*Dcoords%wrk)/REAL(SIZE(Dcoords%wrk, 1)*SIZE(Dcoords%wrk, 2), KIND=dp))
1053 562306 : my_rms_force = SQRT(SUM(forces%wrk*forces%wrk)/REAL(SIZE(forces%wrk, 1)*SIZE(forces%wrk, 2), KIND=dp))
1054 544 : IF (my_max_dr < max_dr) labels(1) = "YES"
1055 544 : IF (my_max_force < max_force) labels(2) = "YES"
1056 544 : IF (my_rms_dr < rms_dr) labels(3) = "YES"
1057 544 : IF (my_rms_force < rms_force) labels(4) = "YES"
1058 586 : IF (ALL(labels == "YES")) converged = .TRUE.
1059 :
1060 : iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "CONVERGENCE_INFO", &
1061 544 : extension=".nebLog")
1062 544 : IF (iw > 0) THEN
1063 : ! Print convergence info
1064 272 : WRITE (iw, FMT='(A,A)') ' **************************************', &
1065 544 : '*****************************************'
1066 : WRITE (iw, FMT='(1X,A,2X,F16.10,5X,"[",F16.10,1X,"]",T76,"(",A,")")') &
1067 272 : 'RMS DISPLACEMENT =', my_rms_dr, rms_dr, labels(3), &
1068 272 : 'MAX DISPLACEMENT =', my_max_dr, max_dr, labels(1), &
1069 272 : 'RMS FORCE =', my_rms_force, rms_force, labels(4), &
1070 544 : 'MAX FORCE =', my_max_force, max_force, labels(2)
1071 272 : WRITE (iw, FMT='(A,A)') ' **************************************', &
1072 544 : '*****************************************'
1073 : END IF
1074 : CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, &
1075 544 : "CONVERGENCE_INFO")
1076 544 : END FUNCTION check_convergence
1077 :
1078 72 : END MODULE neb_utils
|
