LCOV - code coverage report
Current view: top level - src/motion - neb_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 190 199 95.5 %
Date: 2024-03-28 07:31:50 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 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
     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
     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

Generated by: LCOV version 1.15