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 performing a 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 09.2006
19 : !> \par History
20 : !> - Teodoro Laino 10.2008 [tlaino] - University of Zurich
21 : !> Extension to a subspace of collective variables
22 : ! **************************************************************************************************
23 : MODULE neb_methods
24 : USE colvar_utils, ONLY: number_of_colvar
25 : USE cp_external_control, ONLY: external_control
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_type,&
28 : cp_to_string
29 : USE cp_output_handling, ONLY: cp_add_iter_level,&
30 : cp_iterate,&
31 : cp_print_key_finished_output,&
32 : cp_print_key_unit_nr,&
33 : cp_rm_iter_level
34 : USE cp_subsys_types, ONLY: cp_subsys_type
35 : USE f77_interface, ONLY: f_env_add_defaults,&
36 : f_env_rm_defaults,&
37 : f_env_type
38 : USE force_env_types, ONLY: force_env_get
39 : USE global_types, ONLY: global_environment_type
40 : USE header, ONLY: band_header
41 : USE input_constants, ONLY: band_diis_opt,&
42 : band_md_opt,&
43 : do_rep_blocked,&
44 : do_sm
45 : USE input_section_types, ONLY: section_type,&
46 : section_vals_get_subs_vals,&
47 : section_vals_type,&
48 : section_vals_val_get
49 : USE kinds, ONLY: dp
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE neb_io, ONLY: dump_neb_info,&
52 : neb_rep_env_map_info,&
53 : read_neb_section
54 : USE neb_md_utils, ONLY: control_vels_a,&
55 : control_vels_b
56 : USE neb_opt_utils, ONLY: accept_diis_step,&
57 : neb_ls
58 : USE neb_types, ONLY: neb_type,&
59 : neb_var_create,&
60 : neb_var_release,&
61 : neb_var_type
62 : USE neb_utils, ONLY: build_replica_coords,&
63 : check_convergence,&
64 : neb_calc_energy_forces,&
65 : reorient_images,&
66 : reparametrize_images
67 : USE particle_types, ONLY: particle_type
68 : USE physcon, ONLY: massunit
69 : USE replica_methods, ONLY: rep_env_create
70 : USE replica_types, ONLY: rep_env_release,&
71 : replica_env_type
72 : #include "../base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 : PRIVATE
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_methods'
77 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
78 : PUBLIC :: neb
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief Real subroutine for NEB calculations
84 : !> \param input ...
85 : !> \param input_declaration ...
86 : !> \param para_env ...
87 : !> \param globenv ...
88 : !> \author Teodoro Laino 09.2006
89 : !> \note
90 : !> Based on the use of replica_env
91 : ! **************************************************************************************************
92 102 : SUBROUTINE neb(input, input_declaration, para_env, globenv)
93 : TYPE(section_vals_type), POINTER :: input
94 : TYPE(section_type), POINTER :: input_declaration
95 : TYPE(mp_para_env_type), POINTER :: para_env
96 : TYPE(global_environment_type), POINTER :: globenv
97 :
98 : CHARACTER(len=*), PARAMETER :: routineN = 'neb'
99 :
100 : INTEGER :: handle, ierr, iw, iw2, nrep, &
101 : output_unit, prep, proc_dist_type
102 : LOGICAL :: check, row_force
103 : TYPE(cp_logger_type), POINTER :: logger
104 : TYPE(cp_subsys_type), POINTER :: subsys
105 : TYPE(f_env_type), POINTER :: f_env
106 : TYPE(neb_type), POINTER :: neb_env
107 : TYPE(neb_var_type), POINTER :: coords, forces, vels
108 34 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
109 : TYPE(replica_env_type), POINTER :: rep_env
110 : TYPE(section_vals_type), POINTER :: diis_section, force_env_section, &
111 : md_section, motion_section, &
112 : neb_section, print_section
113 :
114 34 : CALL timeset(routineN, handle)
115 34 : NULLIFY (logger, subsys, f_env, rep_env)
116 34 : NULLIFY (forces, coords, vels, neb_env)
117 34 : logger => cp_get_default_logger()
118 34 : CALL cp_add_iter_level(logger%iter_info, "BAND")
119 34 : motion_section => section_vals_get_subs_vals(input, "MOTION")
120 34 : print_section => section_vals_get_subs_vals(motion_section, "PRINT")
121 34 : neb_section => section_vals_get_subs_vals(motion_section, "BAND")
122 : output_unit = cp_print_key_unit_nr(logger, neb_section, "PROGRAM_RUN_INFO", &
123 34 : extension=".nebLog")
124 34 : CALL section_vals_val_get(neb_section, "NPROC_REP", i_val=prep)
125 34 : CALL section_vals_val_get(neb_section, "PROC_DIST_TYPE", i_val=proc_dist_type)
126 34 : row_force = (proc_dist_type == do_rep_blocked)
127 34 : nrep = MAX(1, para_env%num_pe/prep)
128 34 : IF (nrep*prep /= para_env%num_pe .AND. output_unit > 0) THEN
129 : CALL cp_warn(__LOCATION__, &
130 : "Number of totally requested processors ("//TRIM(ADJUSTL(cp_to_string(para_env%num_pe)))//") "// &
131 : "is not compatible with the number of processors requested per replica ("// &
132 : TRIM(ADJUSTL(cp_to_string(prep)))//") and the number of replicas ("// &
133 : TRIM(ADJUSTL(cp_to_string(nrep)))//") . ["// &
134 0 : TRIM(ADJUSTL(cp_to_string(para_env%num_pe - nrep*prep)))//"] processors will be wasted!")
135 : END IF
136 34 : force_env_section => section_vals_get_subs_vals(input, "FORCE_EVAL")
137 : ! Create Replica Environments
138 34 : IF (output_unit > 0) WRITE (output_unit, '(T2,"NEB|",A)') " Replica_env Setup. START"
139 : CALL rep_env_create(rep_env, para_env=para_env, input=input, &
140 34 : input_declaration=input_declaration, nrep=nrep, prep=prep, row_force=row_force)
141 34 : IF (output_unit > 0) WRITE (output_unit, '(T2,"NEB|",A)') " Replica_env Setup. END"
142 34 : IF (ASSOCIATED(rep_env)) THEN
143 34 : CPASSERT(SIZE(rep_env%local_rep_indices) == 1)
144 34 : CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env)
145 34 : CALL force_env_get(f_env%force_env, subsys=subsys)
146 34 : particle_set => subsys%particles%els
147 : ! Read NEB controlling parameters
148 34 : ALLOCATE (neb_env)
149 34 : neb_env%force_env => f_env%force_env
150 34 : neb_env%root_section => input
151 34 : neb_env%force_env_section => force_env_section
152 34 : neb_env%motion_print_section => print_section
153 34 : neb_env%neb_section => neb_section
154 34 : neb_env%nsize_xyz = rep_env%ndim
155 34 : neb_env%nsize_int = number_of_colvar(f_env%force_env)
156 34 : check = (neb_env%nsize_xyz >= neb_env%nsize_int)
157 34 : CPASSERT(check)
158 : ! Check that the used colvar are uniquely determined
159 : check = (number_of_colvar(f_env%force_env) == &
160 34 : number_of_colvar(f_env%force_env, unique=.TRUE.))
161 34 : CPASSERT(check)
162 34 : CALL read_neb_section(neb_env, neb_section)
163 : ! Print BAND header
164 34 : iw2 = cp_print_key_unit_nr(logger, neb_section, "BANNER", extension=".nebLog")
165 34 : CALL band_header(iw2, neb_env%number_of_replica, nrep, prep)
166 34 : CALL cp_print_key_finished_output(iw2, logger, neb_section, "BANNER")
167 : ! Allocate the principal vectors used in the BAND calculation
168 34 : CALL neb_var_create(coords, neb_env, full_allocation=.TRUE.)
169 34 : CALL neb_var_create(forces, neb_env)
170 34 : CALL neb_var_create(vels, neb_env)
171 : ! Collecting the coordinates of the starting replicas of the BAND calculation
172 34 : IF (output_unit > 0) WRITE (output_unit, '(T2,"NEB|",A)') " Building initial set of coordinates. START"
173 : iw = cp_print_key_unit_nr(logger, neb_section, "PROGRAM_RUN_INFO/INITIAL_CONFIGURATION_INFO", &
174 34 : extension=".nebLog")
175 : CALL build_replica_coords(neb_section, particle_set, coords, vels, neb_env, iw, globenv, &
176 34 : rep_env%para_env)
177 : CALL cp_print_key_finished_output(iw, logger, neb_section, &
178 34 : "PROGRAM_RUN_INFO/INITIAL_CONFIGURATION_INFO")
179 34 : IF (output_unit > 0) WRITE (output_unit, '(T2,"NEB|",A)') " Building initial set of coordinates. END"
180 : ! Print some additional info in the replica_env initialization file
181 34 : CALL neb_rep_env_map_info(rep_env, neb_env)
182 : ! Perform NEB optimization
183 50 : SELECT CASE (neb_env%opt_type)
184 : CASE (band_md_opt)
185 16 : neb_env%opt_type_label = "MOLECULAR DYNAMICS"
186 16 : md_section => section_vals_get_subs_vals(neb_section, "OPTIMIZE_BAND%MD")
187 : CALL neb_md(rep_env, neb_env, coords, vels, forces, particle_set, output_unit, &
188 16 : md_section, logger, globenv)
189 : CASE (band_diis_opt)
190 18 : neb_env%opt_type_label = "DIIS"
191 18 : diis_section => section_vals_get_subs_vals(neb_section, "OPTIMIZE_BAND%DIIS")
192 : CALL neb_diis(rep_env, neb_env, coords, vels, forces, particle_set, output_unit, &
193 52 : diis_section, logger, globenv)
194 : END SELECT
195 : ! Release force_eval
196 34 : CALL f_env_rm_defaults(f_env, ierr)
197 : ! Release coords, vels and forces
198 34 : CALL neb_var_release(coords)
199 34 : CALL neb_var_release(forces)
200 34 : CALL neb_var_release(vels)
201 : ! At the end let's destroy the environment of the BAND calculation
202 34 : DEALLOCATE (neb_env)
203 : END IF
204 34 : CALL rep_env_release(rep_env)
205 : CALL cp_print_key_finished_output(output_unit, logger, neb_section, &
206 34 : "PROGRAM_RUN_INFO")
207 34 : CALL cp_rm_iter_level(logger%iter_info, "BAND")
208 34 : CALL timestop(handle)
209 34 : END SUBROUTINE neb
210 :
211 : ! **************************************************************************************************
212 : !> \brief MD type optimization NEB
213 : !> \param rep_env ...
214 : !> \param neb_env ...
215 : !> \param coords ...
216 : !> \param vels ...
217 : !> \param forces ...
218 : !> \param particle_set ...
219 : !> \param output_unit ...
220 : !> \param md_section ...
221 : !> \param logger ...
222 : !> \param globenv ...
223 : !> \author Teodoro Laino 09.2006
224 : ! **************************************************************************************************
225 16 : SUBROUTINE neb_md(rep_env, neb_env, coords, vels, forces, particle_set, output_unit, &
226 : md_section, logger, globenv)
227 : TYPE(replica_env_type), POINTER :: rep_env
228 : TYPE(neb_type), OPTIONAL, POINTER :: neb_env
229 : TYPE(neb_var_type), POINTER :: coords, vels, forces
230 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
231 : INTEGER, INTENT(IN) :: output_unit
232 : TYPE(section_vals_type), POINTER :: md_section
233 : TYPE(cp_logger_type), POINTER :: logger
234 : TYPE(global_environment_type), POINTER :: globenv
235 :
236 : CHARACTER(len=*), PARAMETER :: routineN = 'neb_md'
237 :
238 : INTEGER :: handle, iatom, ic, is, istep, iw, &
239 : max_steps, natom, shell_index
240 : LOGICAL :: converged, should_stop
241 : REAL(KIND=dp) :: dt
242 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: distances, energies
243 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mass
244 : TYPE(neb_var_type), POINTER :: Dcoords
245 : TYPE(section_vals_type), POINTER :: tc_section, vc_section
246 :
247 16 : CALL timeset(routineN, handle)
248 16 : NULLIFY (Dcoords, tc_section, vc_section)
249 16 : CPASSERT(ASSOCIATED(coords))
250 16 : CPASSERT(ASSOCIATED(vels))
251 : ! MD band for string methods type does not make anywa sense. Stop calculation.
252 16 : IF (neb_env%id_type == do_sm) THEN
253 0 : CPWARN("MD band optimization and String Method incompatible.")
254 : END IF
255 : ! Output unit
256 : iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "REPLICA_INFO", &
257 16 : extension=".replicaLog")
258 16 : tc_section => section_vals_get_subs_vals(md_section, "TEMP_CONTROL")
259 16 : vc_section => section_vals_get_subs_vals(md_section, "VEL_CONTROL")
260 16 : CALL section_vals_val_get(md_section, "TIMESTEP", r_val=dt)
261 16 : CALL section_vals_val_get(md_section, "MAX_STEPS", i_val=max_steps)
262 : ! Initial setup for MD
263 16 : CALL neb_var_create(Dcoords, neb_env)
264 64 : ALLOCATE (mass(SIZE(coords%wrk, 1), neb_env%number_of_replica))
265 48 : ALLOCATE (energies(neb_env%number_of_replica))
266 48 : ALLOCATE (distances(neb_env%number_of_replica - 1))
267 : ! Setting up the mass array
268 16 : IF (neb_env%use_colvar) THEN
269 44 : mass(:, :) = 0.5_dp*dt/massunit
270 : ELSE
271 12 : natom = SIZE(particle_set)
272 216 : DO iatom = 1, natom
273 204 : ic = 3*(iatom - 1)
274 204 : shell_index = particle_set(iatom)%shell_index
275 216 : IF (shell_index == 0) THEN
276 4964 : mass(ic + 1:ic + 3, :) = 0.5_dp*dt/particle_set(iatom)%atomic_kind%mass
277 : ELSE
278 0 : is = 3*(natom + shell_index - 1)
279 0 : mass(ic + 1:ic + 3, :) = 0.5_dp*dt/particle_set(iatom)%atomic_kind%shell%mass_core
280 0 : mass(is + 1:is + 3, :) = 0.5_dp*dt/particle_set(iatom)%atomic_kind%shell%mass_shell
281 : END IF
282 : END DO
283 : END IF
284 : ! Initializing forces array
285 : CALL reorient_images(neb_env%rotate_frames, particle_set, coords, vels, &
286 16 : output_unit, distances, neb_env%number_of_replica)
287 90 : neb_env%avg_distance = SQRT(SUM(distances*distances)/REAL(SIZE(distances), KIND=dp))
288 : CALL neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, &
289 16 : particle_set, iw)
290 :
291 : CALL dump_neb_info(neb_env=neb_env, &
292 : coords=coords, &
293 : vels=vels, &
294 : forces=forces, &
295 : particle_set=particle_set, &
296 : logger=logger, &
297 : istep=0, &
298 : energies=energies, &
299 : distances=distances, &
300 16 : output_unit=output_unit)
301 176 : md_opt_loop: DO istep = 1, max_steps
302 160 : CALL cp_iterate(logger%iter_info, iter_nr=istep)
303 : ! Save the optimization step counter
304 160 : neb_env%istep = istep
305 : ! Velocity Verlet (first part)
306 83760 : vels%wrk(:, :) = vels%wrk(:, :) + mass(:, :)*forces%wrk(:, :)
307 : ! Control on velocity - I part [rescale, annealing]
308 : CALL control_vels_a(vels, particle_set, tc_section, vc_section, output_unit, &
309 160 : istep)
310 : ! Coordinate step
311 83760 : Dcoords%wrk(:, :) = dt*vels%wrk(:, :)
312 83760 : coords%wrk(:, :) = coords%wrk(:, :) + Dcoords%wrk(:, :)
313 :
314 : CALL reorient_images(neb_env%rotate_frames, particle_set, coords, vels, &
315 160 : output_unit, distances, neb_env%number_of_replica)
316 900 : neb_env%avg_distance = SQRT(SUM(distances*distances)/REAL(SIZE(distances), KIND=dp))
317 : CALL neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, &
318 160 : particle_set, iw)
319 : ! Check for an external exit command
320 160 : CALL external_control(should_stop, "NEB", globenv=globenv)
321 160 : IF (should_stop) EXIT md_opt_loop
322 : ! Control on velocity - II part [check vels VS forces, Steepest Descent like]
323 160 : CALL control_vels_b(vels, forces, vc_section)
324 : ! Velocity Verlet (second part)
325 83760 : vels%wrk(:, :) = vels%wrk(:, :) + mass(:, :)*forces%wrk(:, :)
326 : ! Dump Infos
327 : CALL dump_neb_info(neb_env=neb_env, &
328 : coords=coords, &
329 : vels=vels, &
330 : forces=forces, &
331 : particle_set=particle_set, &
332 : logger=logger, &
333 : istep=istep, &
334 : energies=energies, &
335 : distances=distances, &
336 160 : output_unit=output_unit)
337 160 : converged = check_convergence(neb_env, Dcoords, forces)
338 336 : IF (converged) EXIT md_opt_loop
339 : END DO md_opt_loop
340 :
341 16 : DEALLOCATE (mass)
342 16 : DEALLOCATE (energies)
343 16 : DEALLOCATE (distances)
344 16 : CALL neb_var_release(Dcoords)
345 : CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, &
346 16 : "REPLICA_INFO")
347 16 : CALL timestop(handle)
348 :
349 32 : END SUBROUTINE neb_md
350 :
351 : ! **************************************************************************************************
352 : !> \brief DIIS type optimization NEB
353 : !> \param rep_env ...
354 : !> \param neb_env ...
355 : !> \param coords ...
356 : !> \param vels ...
357 : !> \param forces ...
358 : !> \param particle_set ...
359 : !> \param output_unit ...
360 : !> \param diis_section ...
361 : !> \param logger ...
362 : !> \param globenv ...
363 : !> \author Teodoro Laino 09.2006
364 : ! **************************************************************************************************
365 18 : SUBROUTINE neb_diis(rep_env, neb_env, coords, vels, forces, particle_set, output_unit, &
366 : diis_section, logger, globenv)
367 : TYPE(replica_env_type), POINTER :: rep_env
368 : TYPE(neb_type), OPTIONAL, POINTER :: neb_env
369 : TYPE(neb_var_type), POINTER :: coords, vels, forces
370 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
371 : INTEGER, INTENT(IN) :: output_unit
372 : TYPE(section_vals_type), POINTER :: diis_section
373 : TYPE(cp_logger_type), POINTER :: logger
374 : TYPE(global_environment_type), POINTER :: globenv
375 :
376 : CHARACTER(len=*), PARAMETER :: routineN = 'neb_diis'
377 :
378 : INTEGER :: handle, istep, iw, iw2, max_sd_steps, &
379 : max_steps, n_diis
380 18 : INTEGER, DIMENSION(:), POINTER :: set_err
381 : LOGICAL :: check_diis, converged, diis_on, do_ls, &
382 : should_stop, skip_ls
383 : REAL(KIND=dp) :: max_stepsize, norm, stepsize, stepsize0
384 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: distances, energies
385 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: crr, err
386 : TYPE(neb_var_type), POINTER :: sline
387 :
388 18 : CALL timeset(routineN, handle)
389 18 : NULLIFY (sline, crr, err)
390 18 : neb_env%opt_type_label = "SD"
391 18 : do_ls = .TRUE.
392 18 : CPASSERT(ASSOCIATED(coords))
393 18 : CPASSERT(ASSOCIATED(vels))
394 18 : CPASSERT(ASSOCIATED(forces))
395 : iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "REPLICA_INFO", &
396 18 : extension=".replicaLog")
397 18 : CALL section_vals_val_get(diis_section, "MAX_STEPS", i_val=max_steps)
398 18 : CALL section_vals_val_get(diis_section, "N_DIIS", i_val=n_diis)
399 18 : CALL section_vals_val_get(diis_section, "STEPSIZE", r_val=stepsize0)
400 18 : CALL section_vals_val_get(diis_section, "MAX_STEPSIZE", r_val=max_stepsize)
401 18 : CALL section_vals_val_get(diis_section, "NO_LS", l_val=skip_ls)
402 18 : CALL section_vals_val_get(diis_section, "MAX_SD_STEPS", i_val=max_sd_steps)
403 18 : CALL section_vals_val_get(diis_section, "CHECK_DIIS", l_val=check_diis)
404 : iw2 = cp_print_key_unit_nr(logger, diis_section, "DIIS_INFO", &
405 18 : extension=".diisLog")
406 : ! Initial setup for DIIS
407 18 : stepsize = stepsize0
408 : ! Allocate type for Line Search direction
409 18 : CALL neb_var_create(sline, neb_env, full_allocation=.TRUE.)
410 : ! Array of error vectors
411 108 : ALLOCATE (err(PRODUCT(coords%size_wrk), n_diis))
412 108 : ALLOCATE (crr(PRODUCT(coords%size_wrk), n_diis))
413 54 : ALLOCATE (set_err(n_diis))
414 54 : ALLOCATE (energies(neb_env%number_of_replica))
415 54 : ALLOCATE (distances(neb_env%number_of_replica - 1))
416 : ! Initializing forces array
417 : CALL reorient_images(neb_env%rotate_frames, particle_set, coords, vels, &
418 18 : output_unit, distances, neb_env%number_of_replica)
419 : CALL reparametrize_images(neb_env%reparametrize_frames, neb_env%spline_order, &
420 18 : neb_env%smoothing, coords%wrk, sline%wrk, distances)
421 122 : neb_env%avg_distance = SQRT(SUM(distances*distances)/REAL(SIZE(distances), KIND=dp))
422 : CALL neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, &
423 18 : particle_set, iw)
424 : ! Dump Infos
425 : CALL dump_neb_info(neb_env=neb_env, &
426 : coords=coords, &
427 : forces=forces, &
428 : particle_set=particle_set, &
429 : logger=logger, &
430 : istep=0, &
431 : energies=energies, &
432 : distances=distances, &
433 : vels=vels, &
434 18 : output_unit=output_unit)
435 : ! If rotation is requested let's apply it at the beginning of the
436 : ! Geometry optimization and then let's disable it
437 18 : neb_env%rotate_frames = .FALSE.
438 : ! Main SD/DIIS loop
439 104 : set_err = -1
440 398 : DO istep = 1, max_steps
441 384 : CALL cp_iterate(logger%iter_info, iter_nr=istep)
442 384 : neb_env%opt_type_label = "SD"
443 : ! Save the optimization step counter
444 384 : neb_env%istep = istep
445 : ! Perform one step of SD with line search
446 520346 : norm = SQRT(SUM(forces%wrk*forces%wrk))
447 384 : IF (norm < EPSILON(0.0_dp)) THEN
448 : ! Let's handle the case in which the band is already fully optimized
449 18 : converged = .TRUE.
450 : EXIT
451 : END IF
452 1040308 : sline%wrk = forces%wrk/norm
453 384 : IF (do_ls .AND. (.NOT. skip_ls)) THEN
454 : CALL neb_ls(stepsize, sline, rep_env, neb_env, coords, energies, forces, &
455 150 : vels, particle_set, iw, output_unit, distances, diis_section, iw2)
456 150 : IF (iw2 > 0) &
457 0 : WRITE (iw2, '(T2,A,T69,F12.6)') "SD| Stepsize in SD after linesearch", &
458 0 : stepsize
459 : ELSE
460 234 : stepsize = MIN(norm*stepsize0, max_stepsize)
461 234 : IF (iw2 > 0) &
462 0 : WRITE (iw2, '(T2,A,T69,F12.6)') "SD| Stepsize in SD no linesearch performed", &
463 0 : stepsize
464 : END IF
465 520346 : sline%wrk = stepsize*sline%wrk
466 : diis_on = accept_diis_step(istep > max_sd_steps, n_diis, err, crr, set_err, sline, coords, &
467 384 : check_diis, iw2)
468 384 : IF (diis_on) THEN
469 146 : neb_env%opt_type_label = "DIIS"
470 : END IF
471 384 : do_ls = .TRUE.
472 1978 : IF (COUNT(set_err == -1) == 1) do_ls = .FALSE.
473 1040308 : coords%wrk = coords%wrk + sline%wrk
474 : ! Compute forces
475 : CALL reorient_images(neb_env%rotate_frames, particle_set, coords, vels, &
476 384 : output_unit, distances, neb_env%number_of_replica)
477 : CALL reparametrize_images(neb_env%reparametrize_frames, neb_env%spline_order, &
478 384 : neb_env%smoothing, coords%wrk, sline%wrk, distances)
479 2462 : neb_env%avg_distance = SQRT(SUM(distances*distances)/REAL(SIZE(distances), KIND=dp))
480 : CALL neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, &
481 384 : particle_set, iw)
482 : ! Check for an external exit command
483 384 : CALL external_control(should_stop, "NEB", globenv=globenv)
484 384 : IF (should_stop) EXIT
485 : ! Dump Infos
486 : CALL dump_neb_info(neb_env=neb_env, &
487 : coords=coords, &
488 : forces=forces, &
489 : particle_set=particle_set, &
490 : logger=logger, &
491 : istep=istep, &
492 : energies=energies, &
493 : distances=distances, &
494 : vels=vels, &
495 384 : output_unit=output_unit)
496 :
497 384 : converged = check_convergence(neb_env, sline, forces)
498 782 : IF (converged) EXIT
499 : END DO
500 18 : DEALLOCATE (energies)
501 18 : DEALLOCATE (distances)
502 18 : DEALLOCATE (err)
503 18 : DEALLOCATE (crr)
504 18 : DEALLOCATE (set_err)
505 18 : CALL neb_var_release(sline)
506 18 : CALL timestop(handle)
507 54 : END SUBROUTINE neb_diis
508 :
509 : END MODULE neb_methods
|
