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