LCOV - code coverage report
Current view: top level - src/motion - neb_opt_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 154 176 87.5 %
Date: 2024-04-20 06:29:22 Functions: 3 3 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 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        1472 :             ALLOCATE (work(np, np))
     105        1472 :             ALLOCATE (wrk_inv(np, np))
     106        1104 :             ALLOCATE (cwrk(np))
     107        1104 :             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        1472 :             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         368 :       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

Generated by: LCOV version 1.15