Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Module with utility to perform MD 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_opt_utils
21 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
22 : cp_to_string
23 : USE input_section_types, ONLY: section_vals_type,&
24 : section_vals_val_get
25 : USE kinds, ONLY: default_string_length,&
26 : dp
27 : USE neb_types, ONLY: neb_type,&
28 : neb_var_type
29 : USE neb_utils, ONLY: neb_calc_energy_forces,&
30 : reorient_images
31 : USE particle_types, ONLY: particle_type
32 : USE replica_types, ONLY: replica_env_type
33 : #include "../base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_opt_utils'
38 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
39 :
40 : PUBLIC :: accept_diis_step, &
41 : neb_ls
42 :
43 : REAL(KIND=dp), DIMENSION(2:10), PRIVATE :: acceptance_factor = &
44 : [0.97_dp, 0.84_dp, 0.71_dp, 0.67_dp, 0.62_dp, 0.56_dp, 0.49_dp, 0.41_dp, 0.0_dp]
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Performs few basic operations for the NEB DIIS optimization
50 : !> \param apply_diis ...
51 : !> \param n_diis ...
52 : !> \param err ...
53 : !> \param crr ...
54 : !> \param set_err ...
55 : !> \param sline ...
56 : !> \param coords ...
57 : !> \param check_diis ...
58 : !> \param iw2 ...
59 : !> \return ...
60 : !> \author Teodoro Laino 10.2006
61 : ! **************************************************************************************************
62 384 : FUNCTION accept_diis_step(apply_diis, n_diis, err, crr, set_err, sline, coords, &
63 : check_diis, iw2) RESULT(accepted)
64 : LOGICAL, INTENT(IN) :: apply_diis
65 : INTEGER, INTENT(IN) :: n_diis
66 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: err, crr
67 : INTEGER, DIMENSION(:), POINTER :: set_err
68 : TYPE(neb_var_type), POINTER :: sline, coords
69 : LOGICAL, INTENT(IN) :: check_diis
70 : INTEGER, INTENT(IN) :: iw2
71 : LOGICAL :: accepted
72 :
73 : CHARACTER(LEN=default_string_length) :: line
74 : INTEGER :: i, iend, ind, indi, indj, info, istart, &
75 : iv, iw, j, jv, k, lwork, np, nv
76 384 : INTEGER, ALLOCATABLE, DIMENSION(:) :: IWORK
77 : LOGICAL :: increase_error
78 384 : REAL(dp), DIMENSION(:, :), POINTER :: work
79 : REAL(KIND=dp) :: eps_svd
80 384 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S, Work_dgesdd
81 384 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: U, VT, wrk, wrk_inv
82 384 : REAL(KIND=dp), DIMENSION(:), POINTER :: awrk, cwrk, ref, step
83 :
84 384 : iw = cp_logger_get_default_io_unit()
85 384 : accepted = .FALSE.
86 : ! find the index with the minimum element of the set_err array
87 1978 : nv = MINLOC(set_err, 1)
88 384 : IF (iw2 > 0) WRITE (iw2, '(A,I3)') "Entering into the DIIS module. Error vector number:: ", nv
89 384 : set_err(nv) = 1
90 384 : eps_svd = 1.0E-10_dp
91 1152 : ALLOCATE (step(sline%size_wrk(1)*sline%size_wrk(2)))
92 1152 : ALLOCATE (ref(sline%size_wrk(1)*sline%size_wrk(2)))
93 518268 : err(:, nv) = RESHAPE(sline%wrk, [sline%size_wrk(1)*sline%size_wrk(2)])
94 518268 : crr(:, nv) = RESHAPE(coords%wrk, [coords%size_wrk(1)*coords%size_wrk(2)])
95 384 : jv = n_diis
96 1658 : IF (ALL(set_err == 1) .AND. apply_diis) THEN
97 186 : IF (iw2 > 0) WRITE (iw2, '(A)') "Applying DIIS equations"
98 : ! Apply DIIS..
99 438 : DO jv = 2, n_diis
100 368 : np = jv + 1
101 368 : IF (iw2 > 0) WRITE (iw2, '(A,I5,A)') "Applying DIIS equations with the last", &
102 0 : jv, " error vectors"
103 1472 : ALLOCATE (wrk(np, np))
104 1104 : ALLOCATE (work(np, np))
105 1104 : ALLOCATE (wrk_inv(np, np))
106 1104 : ALLOCATE (cwrk(np))
107 736 : ALLOCATE (awrk(np))
108 1704 : awrk = 0.0_dp
109 6772 : wrk = 1.0_dp
110 368 : wrk(np, np) = 0.0_dp
111 368 : awrk(np) = 1.0_dp
112 1336 : DO i = 1, jv
113 968 : indi = n_diis - i + 1
114 3202 : DO j = i, jv
115 1866 : indj = n_diis - j + 1
116 2252946 : wrk(i, j) = DOT_PRODUCT(err(:, indi), err(:, indj))
117 2834 : wrk(j, i) = wrk(i, j)
118 : END DO
119 : END DO
120 368 : IF (iw2 > 0) THEN
121 0 : line = "DIIS Matrix"//cp_to_string(np)//"x"//cp_to_string(np)//"."
122 0 : WRITE (iw2, '(A)') TRIM(line)
123 0 : WRITE (iw2, '('//cp_to_string(np)//'F12.6)') wrk
124 : END IF
125 : ! Inverte the DIIS Matrix
126 6772 : work = TRANSPOSE(wrk)
127 : ! Workspace query
128 1104 : ALLOCATE (iwork(8*np))
129 1104 : ALLOCATE (S(np))
130 1472 : ALLOCATE (U(np, np))
131 1104 : ALLOCATE (VT(np, np))
132 368 : ALLOCATE (work_dgesdd(1))
133 368 : lwork = -1
134 368 : CALL DGESDD('S', np, np, work, np, S, U, np, vt, np, work_dgesdd, lwork, iwork, info)
135 368 : lwork = INT(work_dgesdd(1))
136 368 : DEALLOCATE (work_dgesdd)
137 1104 : ALLOCATE (work_dgesdd(lwork))
138 368 : CALL DGESDD('S', np, np, work, np, S, U, np, vt, np, work_dgesdd, lwork, iwork, info)
139 : ! Construct the inverse
140 1704 : DO k = 1, np
141 : ! Invert SV
142 1336 : IF (S(k) < eps_svd) THEN
143 0 : S(k) = 0.0_dp
144 : ELSE
145 1336 : S(k) = 1.0_dp/S(k)
146 : END IF
147 6772 : VT(k, :) = VT(k, :)*S(k)
148 : END DO
149 368 : CALL DGEMM('T', 'T', np, np, np, 1.0_dp, VT, np, U, np, 0.0_dp, wrk_inv, np)
150 368 : DEALLOCATE (iwork)
151 368 : DEALLOCATE (S)
152 368 : DEALLOCATE (U)
153 368 : DEALLOCATE (VT)
154 368 : DEALLOCATE (work_dgesdd)
155 9444 : cwrk = MATMUL(wrk_inv, awrk)
156 : ! Check the DIIS solution
157 469760 : step = 0.0_dp
158 368 : ind = 0
159 1336 : DO i = n_diis, n_diis - jv + 1, -1
160 968 : ind = ind + 1
161 2413420 : step = step + (crr(:, i) + err(:, i))*cwrk(ind)
162 : END DO
163 939152 : step = step - crr(:, n_diis)
164 939152 : ref = err(:, n_diis)
165 : increase_error = check_diis_solution(jv, cwrk, step, ref, &
166 368 : iw2, check_diis)
167 : ! possibly enlarge the error space
168 368 : IF (increase_error) THEN
169 252 : accepted = .TRUE.
170 283380 : sline%wrk = RESHAPE(step, [sline%size_wrk(1), sline%size_wrk(2)])
171 : ELSE
172 116 : DEALLOCATE (awrk)
173 116 : DEALLOCATE (cwrk)
174 116 : DEALLOCATE (wrk)
175 116 : DEALLOCATE (work)
176 116 : DEALLOCATE (wrk_inv)
177 : EXIT
178 : END IF
179 252 : DEALLOCATE (awrk)
180 252 : DEALLOCATE (cwrk)
181 252 : DEALLOCATE (wrk)
182 252 : DEALLOCATE (work)
183 322 : DEALLOCATE (wrk_inv)
184 : END DO
185 186 : IF (iw2 > 0) THEN
186 0 : line = "Exiting DIIS accepting"//cp_to_string(MIN(n_diis, jv))//" errors."
187 0 : WRITE (iw2, '(A)') TRIM(line)
188 : END IF
189 : END IF
190 1658 : IF (ALL(set_err == 1)) THEN
191 : ! always delete the last error vector from the history vectors
192 : ! move error vectors and the set_err in order to have free space
193 : ! at the end of the err array
194 198 : istart = MAX(2, n_diis - jv + 2)
195 198 : iend = n_diis
196 198 : indi = 0
197 626 : DO iv = istart, iend
198 428 : indi = indi + 1
199 509420 : err(:, indi) = err(:, iv)
200 509420 : crr(:, indi) = crr(:, iv)
201 626 : set_err(indi) = 1
202 : END DO
203 530 : DO iv = indi + 1, iend
204 530 : set_err(iv) = -1
205 : END DO
206 : END IF
207 384 : DEALLOCATE (step)
208 384 : DEALLOCATE (ref)
209 :
210 768 : END FUNCTION accept_diis_step
211 :
212 : ! **************************************************************************************************
213 : !> \brief Check conditions for the acceptance of the DIIS step
214 : !> \param nv ...
215 : !> \param cwrk ...
216 : !> \param step ...
217 : !> \param ref ...
218 : !> \param output_unit ...
219 : !> \param check_diis ...
220 : !> \return ...
221 : !> \author Teodoro Laino 10.2006
222 : ! **************************************************************************************************
223 368 : FUNCTION check_diis_solution(nv, cwrk, step, ref, output_unit, check_diis) &
224 : RESULT(accepted)
225 : INTEGER, INTENT(IN) :: nv
226 : REAL(KIND=dp), DIMENSION(:), POINTER :: cwrk, step, ref
227 : INTEGER, INTENT(IN) :: output_unit
228 : LOGICAL, INTENT(IN) :: check_diis
229 : LOGICAL :: accepted
230 :
231 : REAL(KIND=dp) :: costh, norm1, norm2
232 368 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp
233 :
234 368 : accepted = .TRUE.
235 1104 : ALLOCATE (tmp(SIZE(step)))
236 : IF (accepted) THEN
237 : ! (a) The direction of the DIIS step, can be compared to the reference step.
238 : ! if the angle is grater than a specified value, the DIIS step is not
239 : ! acceptable.
240 469760 : norm1 = SQRT(DOT_PRODUCT(ref, ref))
241 469760 : norm2 = SQRT(DOT_PRODUCT(step, step))
242 469760 : costh = DOT_PRODUCT(ref, step)/(norm1*norm2)
243 368 : IF (check_diis) THEN
244 320 : IF (costh < acceptance_factor(MIN(10, nv))) accepted = .FALSE.
245 : ELSE
246 48 : IF (costh <= 0.0_dp) accepted = .FALSE.
247 : END IF
248 368 : IF (output_unit > 0 .AND. (.NOT. accepted)) THEN
249 : WRITE (output_unit, '(T2,"DIIS|",A)') &
250 0 : "The direction of the DIIS step, can be compared to the reference step.", &
251 0 : "If the angle is grater than a specified value, the DIIS step is not", &
252 0 : "acceptable. Value exceeded. Reset DIIS!"
253 : WRITE (output_unit, '(T2,"DIIS|",A,F6.3,A,F6.3,A)') &
254 0 : "Present Cosine <", costh, "> compared with the optimal value <", &
255 0 : acceptance_factor(MIN(10, nv)), "> ."
256 : END IF
257 : END IF
258 368 : IF (accepted .AND. check_diis) THEN
259 : ! (b) The length of the DIIS step is limited to be no more than 10 times
260 : ! the reference step
261 212 : IF (norm1 > norm2*10.0_dp) accepted = .FALSE.
262 212 : IF (output_unit > 0 .AND. (.NOT. accepted)) THEN
263 : WRITE (output_unit, '(T2,"DIIS|",A)') &
264 0 : "The length of the DIIS step is limited to be no more than 10 times", &
265 0 : "the reference step. Value exceeded. Reset DIIS!"
266 : END IF
267 : END IF
268 252 : IF (accepted .AND. check_diis) THEN
269 : ! (d) If the DIIS matrix is nearly singular, the norm of the DIIS step
270 : ! vector becomes small and cwrk/norm1 becomes large, signaling a
271 : ! numerical stability problems. If the magnitude of cwrk/norm1
272 : ! exceeds 10^8 then the step size is assumed to be unacceptable
273 730 : IF (ANY(ABS(cwrk(1:nv)/norm1) > 10**8_dp)) accepted = .FALSE.
274 212 : IF (output_unit > 0 .AND. (.NOT. accepted)) THEN
275 : WRITE (output_unit, '(T2,"DIIS|",A)') &
276 0 : "If the DIIS matrix is nearly singular, the norm of the DIIS step", &
277 0 : "vector becomes small and Coeff/E_norm becomes large, signaling a", &
278 0 : "numerical stability problems. IF the magnitude of Coeff/E_norm", &
279 0 : "exceeds 10^8 THEN the step size is assumed to be unacceptable.", &
280 0 : "Value exceeded. Reset DIIS!"
281 : END IF
282 : END IF
283 368 : DEALLOCATE (tmp)
284 368 : END FUNCTION check_diis_solution
285 :
286 : ! **************************************************************************************************
287 : !> \brief Perform a line minimization search in optimizing a NEB with DIIS
288 : !> \param stepsize ...
289 : !> \param sline ...
290 : !> \param rep_env ...
291 : !> \param neb_env ...
292 : !> \param coords ...
293 : !> \param energies ...
294 : !> \param forces ...
295 : !> \param vels ...
296 : !> \param particle_set ...
297 : !> \param iw ...
298 : !> \param output_unit ...
299 : !> \param distances ...
300 : !> \param diis_section ...
301 : !> \param iw2 ...
302 : !> \author Teodoro Laino 10.2006
303 : ! **************************************************************************************************
304 450 : SUBROUTINE neb_ls(stepsize, sline, rep_env, neb_env, coords, energies, forces, &
305 150 : vels, particle_set, iw, output_unit, distances, diis_section, iw2)
306 : REAL(KIND=dp), INTENT(INOUT) :: stepsize
307 : TYPE(neb_var_type), POINTER :: sline
308 : TYPE(replica_env_type), POINTER :: rep_env
309 : TYPE(neb_type), OPTIONAL, POINTER :: neb_env
310 : TYPE(neb_var_type), POINTER :: coords
311 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: energies
312 : TYPE(neb_var_type), POINTER :: forces, vels
313 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
314 : INTEGER, INTENT(IN) :: iw, output_unit
315 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: distances
316 : TYPE(section_vals_type), POINTER :: diis_section
317 : INTEGER, INTENT(IN) :: iw2
318 :
319 : INTEGER :: i, np
320 : REAL(KIND=dp) :: a, b, max_stepsize, xa, xb, xc_cray, ya, &
321 : yb
322 150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Icoord
323 :
324 : ! replaced xc by xc_cray to work around yet another bug in pgf90 on CRAY xt3
325 :
326 600 : ALLOCATE (Icoord(coords%size_wrk(1), coords%size_wrk(2)))
327 150 : CALL section_vals_val_get(diis_section, "NP_LS", i_val=np)
328 150 : CALL section_vals_val_get(diis_section, "MAX_STEPSIZE", r_val=max_stepsize)
329 247768 : Icoord(:, :) = coords%wrk
330 150 : xa = 0.0_dp
331 247768 : ya = SUM(sline%wrk*forces%wrk)
332 150 : xb = xa + MIN(ya*stepsize, max_stepsize)
333 150 : xc_cray = xb
334 150 : i = 1
335 462 : DO WHILE (i <= np - 1)
336 330 : i = i + 1
337 961582 : coords%wrk = Icoord + xb*sline%wrk
338 : CALL reorient_images(neb_env%rotate_frames, particle_set, coords, vels, &
339 330 : output_unit, distances, neb_env%number_of_replica)
340 1466 : neb_env%avg_distance = SQRT(SUM(distances*distances)/REAL(SIZE(distances), KIND=dp))
341 : CALL neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, &
342 330 : particle_set, iw)
343 480956 : yb = SUM(sline%wrk*forces%wrk)
344 330 : a = (ya - yb)/(2.0_dp*(xa - xb))
345 330 : b = ya - 2.0_dp*a*xa
346 330 : xc_cray = -b/(2.0_dp*a)
347 330 : IF (xc_cray > max_stepsize) THEN
348 18 : IF (iw2 > 0) WRITE (iw2, '(T2,2(A,F6.3),A)') &
349 0 : "LS| Predicted stepsize (", xc_cray, ") greater than allowed stepsize (", &
350 0 : max_stepsize, "). Reset!"
351 18 : xc_cray = max_stepsize
352 18 : EXIT
353 : END IF
354 : ! No Extrapolation .. only interpolation
355 312 : IF ((xc_cray <= MIN(xa, xb) .OR. xc_cray >= MAX(xa, xb)) .AND. (ABS(xa - xb) > 1.0E-5_dp)) THEN
356 180 : IF (iw2 > 0) WRITE (iw2, '(T2,2(A,I5),A)') &
357 0 : "LS| Increasing the number of point from ", np, " to ", np + 1, "."
358 180 : np = np + 1
359 : END IF
360 : !
361 312 : IF (ABS(yb) < ABS(ya)) THEN
362 290 : ya = yb
363 290 : xa = xb
364 : END IF
365 132 : xb = xc_cray
366 : END DO
367 150 : stepsize = xc_cray
368 247768 : coords%wrk = Icoord
369 150 : DEALLOCATE (Icoord)
370 150 : END SUBROUTINE neb_ls
371 :
372 368 : END MODULE neb_opt_utils
|
