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