LCOV - code coverage report
Current view: top level - src/motion - bfgs_optimizer.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 515 548 94.0 %
Date: 2024-04-20 06:29:22 Functions: 12 12 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 Routines for Geometry optimization using BFGS algorithm
      10             : !> \par History
      11             : !>      Module modified by Pierre-André Cazade [pcazade] 01.2020 - University of Limerick.
      12             : !>      Modifications enable Space Group Symmetry.
      13             : ! **************************************************************************************************
      14             : MODULE bfgs_optimizer
      15             : 
      16             :    USE atomic_kind_list_types, ONLY: atomic_kind_list_type
      17             :    USE atomic_kind_types, ONLY: get_atomic_kind, &
      18             :                                 get_atomic_kind_set
      19             :    USE cell_types, ONLY: cell_type, &
      20             :                          pbc
      21             :    USE constraint_fxd, ONLY: fix_atom_control
      22             :    USE cp_blacs_env, ONLY: cp_blacs_env_create, &
      23             :                            cp_blacs_env_release, &
      24             :                            cp_blacs_env_type
      25             :    USE cp_external_control, ONLY: external_control
      26             :    USE cp_files, ONLY: close_file, &
      27             :                        open_file
      28             :    USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale, &
      29             :                                  cp_fm_transpose
      30             :    USE cp_fm_diag, ONLY: choose_eigv_solver
      31             :    USE cp_fm_struct, ONLY: cp_fm_struct_create, &
      32             :                            cp_fm_struct_release, &
      33             :                            cp_fm_struct_type
      34             :    USE cp_fm_types, ONLY: &
      35             :       cp_fm_get_info, &
      36             :       cp_fm_read_unformatted, &
      37             :       cp_fm_set_all, &
      38             :       cp_fm_to_fm, &
      39             :       cp_fm_type, &
      40             :       cp_fm_write_unformatted, cp_fm_create, cp_fm_release
      41             :    USE parallel_gemm_api, ONLY: parallel_gemm
      42             :    USE cp_log_handling, ONLY: cp_get_default_logger, &
      43             :                               cp_logger_type, &
      44             :                               cp_to_string
      45             :    USE cp_output_handling, ONLY: cp_iterate, &
      46             :                                  cp_p_file, &
      47             :                                  cp_print_key_finished_output, &
      48             :                                  cp_print_key_should_output, &
      49             :                                  cp_print_key_unit_nr
      50             :    USE message_passing, ONLY: mp_para_env_type
      51             :    USE cp_subsys_types, ONLY: cp_subsys_get, &
      52             :                               cp_subsys_type
      53             :    USE force_env_types, ONLY: force_env_get, &
      54             :                               force_env_type
      55             :    USE global_types, ONLY: global_environment_type
      56             :    USE gopt_f_methods, ONLY: gopt_f_ii, &
      57             :                              gopt_f_io, &
      58             :                              gopt_f_io_finalize, &
      59             :                              gopt_f_io_init, &
      60             :                              print_geo_opt_header, &
      61             :                              print_geo_opt_nc
      62             :    USE gopt_f_types, ONLY: gopt_f_type
      63             :    USE gopt_param_types, ONLY: gopt_param_type
      64             :    USE input_constants, ONLY: default_cell_method_id, &
      65             :                               default_ts_method_id
      66             :    USE input_section_types, ONLY: section_vals_get_subs_vals, &
      67             :                                   section_vals_type, &
      68             :                                   section_vals_val_get, &
      69             :                                   section_vals_val_set
      70             :    USE kinds, ONLY: default_path_length, &
      71             :                     dp
      72             :    USE machine, ONLY: m_flush, &
      73             :                       m_walltime
      74             :    USE particle_list_types, ONLY: particle_list_type
      75             :    USE space_groups, ONLY: identify_space_group, &
      76             :                            print_spgr, &
      77             :                            spgr_apply_rotations_coord, &
      78             :                            spgr_apply_rotations_force
      79             :    USE space_groups_types, ONLY: spgr_type
      80             : #include "../base/base_uses.f90"
      81             : 
      82             :    IMPLICIT NONE
      83             :    PRIVATE
      84             : 
      85             :    #:include "gopt_f77_methods.fypp"
      86             : 
      87             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bfgs_optimizer'
      88             :    LOGICAL, PARAMETER                   :: debug_this_module = .TRUE.
      89             : 
      90             :    PUBLIC :: geoopt_bfgs
      91             : 
      92             : CONTAINS
      93             : 
      94             : ! **************************************************************************************************
      95             : !> \brief Main driver for BFGS geometry optimizations
      96             : !> \param force_env ...
      97             : !> \param gopt_param ...
      98             : !> \param globenv ...
      99             : !> \param geo_section ...
     100             : !> \param gopt_env ...
     101             : !> \param x0 ...
     102             : !> \par History
     103             : !>      01.2020 modified to perform Space Group Symmetry [pcazade]
     104             : ! **************************************************************************************************
     105        1260 :    RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, gopt_env, x0)
     106             : 
     107             :       TYPE(force_env_type), POINTER                      :: force_env
     108             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     109             :       TYPE(global_environment_type), POINTER             :: globenv
     110             :       TYPE(section_vals_type), POINTER                   :: geo_section
     111             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     112             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
     113             : 
     114             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'geoopt_bfgs'
     115             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     116             : 
     117             :       CHARACTER(LEN=5)                                   :: wildcard
     118             :       CHARACTER(LEN=default_path_length)                 :: hes_filename
     119             :       INTEGER                                            :: handle, hesunit_read, indf, info, &
     120             :                                                             iter_nr, its, maxiter, ndf, nfree, &
     121             :                                                             output_unit
     122             :       LOGICAL                                            :: conv, hesrest, ionode, shell_present, &
     123             :                                                             should_stop, use_mod_hes, use_rfo
     124             :       REAL(KIND=dp)                                      :: ediff, emin, eold, etot, pred, rad, rat, &
     125             :                                                             step, t_diff, t_now, t_old
     126        1260 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: dg, dr, dx, eigval, gold, work, xold
     127        1260 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: g
     128             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     129             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     130             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_hes
     131             :       TYPE(cp_fm_type)                          :: eigvec_mat, hess_mat, hess_tmp
     132             :       TYPE(cp_logger_type), POINTER                      :: logger
     133             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     134             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     135             :       TYPE(section_vals_type), POINTER                   :: print_key, root_section
     136             :       TYPE(spgr_type), POINTER                           :: spgr
     137             : 
     138        1260 :       NULLIFY (logger, g, blacs_env, spgr)
     139        2520 :       logger => cp_get_default_logger()
     140        1260 :       para_env => force_env%para_env
     141        1260 :       root_section => force_env%root_section
     142        1260 :       spgr => gopt_env%spgr
     143        1260 :       t_old = m_walltime()
     144             : 
     145        1260 :       CALL timeset(routineN, handle)
     146        1260 :       CALL section_vals_val_get(geo_section, "BFGS%TRUST_RADIUS", r_val=rad)
     147        1260 :       print_key => section_vals_get_subs_vals(geo_section, "BFGS%RESTART")
     148        1260 :       ionode = para_env%is_source()
     149        1260 :       maxiter = gopt_param%max_iter
     150        1260 :       conv = .FALSE.
     151        1260 :       rat = 0.0_dp
     152        1260 :       wildcard = " BFGS"
     153             : 
     154             :       ! Stop if not yet implemented
     155        1260 :       SELECT CASE (gopt_env%type_id)
     156             :       CASE (default_ts_method_id)
     157        1260 :          CPABORT("BFGS method not yet working with DIMER")
     158             :       END SELECT
     159             : 
     160        1260 :       CALL section_vals_val_get(geo_section, "BFGS%USE_RAT_FUN_OPT", l_val=use_rfo)
     161        1260 :       CALL section_vals_val_get(geo_section, "BFGS%USE_MODEL_HESSIAN", l_val=use_mod_hes)
     162        1260 :       CALL section_vals_val_get(geo_section, "BFGS%RESTART_HESSIAN", l_val=hesrest)
     163             :       output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", &
     164        1260 :                                          extension=".geoLog")
     165        1260 :       IF (output_unit > 0) THEN
     166         630 :          IF (use_rfo) THEN
     167             :             WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") &
     168          35 :                "BFGS| Use rational function optimization for step estimation: ", "YES"
     169             :          ELSE
     170             :             WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") &
     171         595 :                "BFGS| Use rational function optimization for step estimation: ", " NO"
     172             :          END IF
     173         630 :          IF (use_mod_hes) THEN
     174             :             WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
     175         526 :                "BFGS| Use model Hessian for initial guess: ", "YES"
     176             :          ELSE
     177             :             WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
     178         104 :                "BFGS| Use model Hessian for initial guess: ", " NO"
     179             :          END IF
     180         630 :          IF (hesrest) THEN
     181             :             WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
     182           1 :                "BFGS| Restart Hessian: ", "YES"
     183             :          ELSE
     184             :             WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
     185         629 :                "BFGS| Restart Hessian: ", " NO"
     186             :          END IF
     187             :          WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.3)") &
     188         630 :             "BFGS| Trust radius: ", rad
     189             :       END IF
     190             : 
     191        1260 :       ndf = SIZE(x0)
     192        1260 :       nfree = gopt_env%nfree
     193        1260 :       IF (ndf > 3000) &
     194             :          CALL cp_warn(__LOCATION__, &
     195             :                       "The dimension of the Hessian matrix ("// &
     196             :                       TRIM(ADJUSTL(cp_to_string(ndf)))//") is greater than 3000. "// &
     197             :                       "The diagonalisation of the full Hessian  matrix needed for BFGS "// &
     198             :                       "is computationally expensive. You should consider to use the linear "// &
     199           0 :                       "scaling variant L-BFGS instead.")
     200             : 
     201             :       ! Initialize hessian (hes = unitary matrix or model hessian )
     202             :       CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
     203        1260 :                                globenv%blacs_repeatable)
     204             :       CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
     205        1260 :                                nrow_global=ndf, ncol_global=ndf)
     206        1260 :       CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
     207        1260 :       CALL cp_fm_create(hess_tmp, fm_struct_hes, name="hess_tmp")
     208        1260 :       CALL cp_fm_create(eigvec_mat, fm_struct_hes, name="eigvec_mat")
     209        3780 :       ALLOCATE (eigval(ndf))
     210       89700 :       eigval(:) = 0.0_dp
     211             : 
     212        1260 :       CALL force_env_get(force_env=force_env, subsys=subsys)
     213        1260 :       CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
     214        1260 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
     215        1260 :       IF (use_mod_hes) THEN
     216        1052 :          IF (shell_present) THEN
     217             :             CALL cp_warn(__LOCATION__, &
     218             :                          "No model Hessian is available for core-shell models. "// &
     219           4 :                          "A unit matrix is used as the initial Hessian.")
     220           4 :             use_mod_hes = .FALSE.
     221             :          END IF
     222        1052 :          IF (gopt_env%type_id == default_cell_method_id) THEN
     223             :             CALL cp_warn(__LOCATION__, &
     224             :                          "No model Hessian is available for cell optimizations. "// &
     225           0 :                          "A unit matrix is used as the initial Hessian.")
     226           0 :             use_mod_hes = .FALSE.
     227             :          END IF
     228             :       END IF
     229             : 
     230        1260 :       IF (use_mod_hes) THEN
     231        1048 :          CALL cp_fm_set_all(hess_mat, alpha=zero, beta=0.00_dp)
     232        1048 :          CALL construct_initial_hess(gopt_env%force_env, hess_mat)
     233        1048 :          CALL cp_fm_to_fm(hess_mat, hess_tmp)
     234        1048 :          CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval, info=info)
     235             :          ! In rare cases the diagonalization of hess_mat fails (bug in scalapack?)
     236        1048 :          IF (info /= 0) THEN
     237           0 :             CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
     238           0 :             IF (output_unit > 0) WRITE (output_unit, *) &
     239           0 :                "BFGS: Matrix diagonalization failed, using unity as model Hessian."
     240             :          ELSE
     241       55822 :             DO its = 1, SIZE(eigval)
     242       55822 :                IF (eigval(its) .LT. 0.1_dp) eigval(its) = 0.1_dp
     243             :             END DO
     244        1048 :             CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
     245        1048 :             CALL cp_fm_column_scale(eigvec_mat, eigval)
     246        1048 :             CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, hess_mat)
     247             :          END IF
     248             :       ELSE
     249         212 :          CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
     250             :       END IF
     251             : 
     252        3780 :       ALLOCATE (xold(ndf))
     253       89700 :       xold(:) = x0(:)
     254             : 
     255        3780 :       ALLOCATE (g(ndf))
     256       89700 :       g(:) = 0.0_dp
     257             : 
     258        3780 :       ALLOCATE (gold(ndf))
     259       89700 :       gold(:) = 0.0_dp
     260             : 
     261        3780 :       ALLOCATE (dx(ndf))
     262       89700 :       dx(:) = 0.0_dp
     263             : 
     264        3780 :       ALLOCATE (dg(ndf))
     265       89700 :       dg(:) = 0.0_dp
     266             : 
     267        3780 :       ALLOCATE (work(ndf))
     268       89700 :       work(:) = 0.0_dp
     269             : 
     270        3780 :       ALLOCATE (dr(ndf))
     271       89700 :       dr(:) = 0.0_dp
     272             : 
     273             :       ! find space_group
     274        1260 :       CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
     275        1260 :       IF (spgr%keep_space_group) THEN
     276           4 :          CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
     277           4 :          CALL spgr_apply_rotations_coord(spgr, x0)
     278           4 :          CALL print_spgr(spgr)
     279             :       END IF
     280             : 
     281             :       ! Geometry optimization starts now
     282        1260 :       CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
     283        1260 :       CALL print_geo_opt_header(gopt_env, output_unit, wildcard)
     284             : 
     285             :       ! Calculate Energy & Gradients
     286             :       CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, &
     287        1260 :                       .FALSE., gopt_env%force_env%para_env)
     288             : 
     289             :       ! Symmetrize coordinates and forces
     290        1260 :       IF (spgr%keep_space_group) THEN
     291           4 :          CALL spgr_apply_rotations_coord(spgr, x0)
     292           4 :          CALL spgr_apply_rotations_force(spgr, g)
     293             :       END IF
     294             : 
     295             :       ! Print info at time 0
     296        1260 :       emin = etot
     297        1260 :       t_now = m_walltime()
     298        1260 :       t_diff = t_now - t_old
     299        1260 :       t_old = t_now
     300        1260 :       CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff)
     301        7367 :       DO its = iter_nr + 1, maxiter
     302        7357 :          CALL cp_iterate(logger%iter_info, last=(its == maxiter))
     303        7357 :          CALL section_vals_val_set(geo_section, "STEP_START_VAL", i_val=its)
     304        7357 :          CALL gopt_f_ii(its, output_unit)
     305             : 
     306             :          ! Hessian update/restarting
     307        7357 :          IF (((its - iter_nr) == 1) .AND. hesrest) THEN
     308           2 :             IF (ionode) THEN
     309           1 :                CALL section_vals_val_get(geo_section, "BFGS%RESTART_FILE_NAME", c_val=hes_filename)
     310             :                CALL open_file(file_name=hes_filename, file_status="OLD", &
     311           1 :                               file_form="UNFORMATTED", file_action="READ", unit_number=hesunit_read)
     312             :             END IF
     313           2 :             CALL cp_fm_read_unformatted(hess_mat, hesunit_read)
     314           2 :             IF (ionode) CALL close_file(unit_number=hesunit_read)
     315             :          ELSE
     316        7355 :             IF ((its - iter_nr) > 1) THEN
     317             :                ! Symmetrize old coordinates and old forces
     318        6107 :                IF (spgr%keep_space_group) THEN
     319           0 :                   CALL spgr_apply_rotations_coord(spgr, xold)
     320           0 :                   CALL spgr_apply_rotations_force(spgr, gold)
     321             :                END IF
     322             : 
     323      606722 :                DO indf = 1, ndf
     324      600615 :                   dx(indf) = x0(indf) - xold(indf)
     325      606722 :                   dg(indf) = g(indf) - gold(indf)
     326             :                END DO
     327             : 
     328        6107 :                CALL bfgs(ndf, dx, dg, hess_mat, work, para_env)
     329             : 
     330             :                ! Symmetrize coordinates and forces change
     331        6107 :                IF (spgr%keep_space_group) THEN
     332           0 :                   CALL spgr_apply_rotations_force(spgr, dx)
     333           0 :                   CALL spgr_apply_rotations_force(spgr, dg)
     334             :                END IF
     335             : 
     336             :                !Possibly dump the Hessian file
     337        6107 :                IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     338        3763 :                   CALL write_bfgs_hessian(geo_section, hess_mat, logger)
     339             :                END IF
     340             :             END IF
     341             :          END IF
     342             : 
     343             :          ! Symmetrize coordinates and forces
     344        7357 :          IF (spgr%keep_space_group) THEN
     345           4 :             CALL spgr_apply_rotations_coord(spgr, x0)
     346           4 :             CALL spgr_apply_rotations_force(spgr, g)
     347             :          END IF
     348             : 
     349             :          ! Setting the present positions & gradients as old
     350      696268 :          xold(:) = x0
     351      696268 :          gold(:) = g
     352             : 
     353             :          ! Copying hessian hes to (ndf x ndf) matrix hes_mat for diagonalization
     354        7357 :          CALL cp_fm_to_fm(hess_mat, hess_tmp)
     355             : 
     356        7357 :          CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval, info=info)
     357             : 
     358             :          ! In rare cases the diagonalization of hess_mat fails (bug in scalapack?)
     359        7357 :          IF (info /= 0) THEN
     360           0 :             IF (output_unit > 0) WRITE (output_unit, *) &
     361           0 :                "BFGS: Matrix diagonalization failed, resetting Hessian to unity."
     362           0 :             CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
     363           0 :             CALL cp_fm_to_fm(hess_mat, hess_tmp)
     364           0 :             CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval)
     365             :          END IF
     366             : 
     367        7357 :          IF (use_rfo) THEN
     368        1712 :             CALL set_hes_eig(ndf, eigval, work)
     369      131300 :             dx(:) = eigval
     370        1712 :             CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
     371             :          END IF
     372        7357 :          CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
     373             : 
     374             :          ! Symmetrize dr
     375        7357 :          IF (spgr%keep_space_group) THEN
     376           4 :             CALL spgr_apply_rotations_force(spgr, dr)
     377             :          END IF
     378             : 
     379        7357 :          CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
     380             : 
     381             :          ! Update the atomic positions
     382      696268 :          x0 = x0 + dr
     383             : 
     384             :          ! Symmetrize coordinates
     385        7357 :          IF (spgr%keep_space_group) THEN
     386           4 :             CALL spgr_apply_rotations_coord(spgr, x0)
     387             :          END IF
     388             : 
     389        7357 :          CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
     390        7357 :          eold = etot
     391             : 
     392             :          ! Energy & Gradients at new step
     393             :          CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, &
     394        7357 :                          .FALSE., gopt_env%force_env%para_env)
     395             : 
     396        7357 :          ediff = etot - eold
     397             : 
     398             :          ! Symmetrize forces
     399        7357 :          IF (spgr%keep_space_group) THEN
     400           4 :             CALL spgr_apply_rotations_force(spgr, g)
     401             :          END IF
     402             : 
     403             :          ! check for an external exit command
     404        7357 :          CALL external_control(should_stop, "GEO", globenv=globenv)
     405        7357 :          IF (should_stop) EXIT
     406             : 
     407             :          ! Some IO and Convergence check
     408        7357 :          t_now = m_walltime()
     409        7357 :          t_diff = t_now - t_old
     410        7357 :          t_old = t_now
     411             :          CALL gopt_f_io(gopt_env, force_env, root_section, its, etot, output_unit, &
     412             :                         eold, emin, wildcard, gopt_param, ndf, dr, g, conv, pred, rat, &
     413        7357 :                         step, rad, used_time=t_diff)
     414             : 
     415        7357 :          IF (conv .OR. (its == maxiter)) EXIT
     416        6107 :          IF (etot < emin) emin = etot
     417       22081 :          IF (use_rfo) CALL update_trust_rad(rat, rad, step, ediff)
     418             :       END DO
     419             : 
     420        1260 :       IF (its == maxiter .AND. (.NOT. conv)) THEN
     421         599 :          CALL print_geo_opt_nc(gopt_env, output_unit)
     422             :       END IF
     423             : 
     424             :       ! Write final  information, if converged
     425        1260 :       CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
     426        1260 :       CALL write_bfgs_hessian(geo_section, hess_mat, logger)
     427             :       CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
     428        1260 :                               gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
     429             : 
     430        1260 :       CALL cp_fm_struct_release(fm_struct_hes)
     431        1260 :       CALL cp_fm_release(hess_mat)
     432        1260 :       CALL cp_fm_release(eigvec_mat)
     433        1260 :       CALL cp_fm_release(hess_tmp)
     434             : 
     435        1260 :       CALL cp_blacs_env_release(blacs_env)
     436        1260 :       DEALLOCATE (xold)
     437        1260 :       DEALLOCATE (g)
     438        1260 :       DEALLOCATE (gold)
     439        1260 :       DEALLOCATE (dx)
     440        1260 :       DEALLOCATE (dg)
     441        1260 :       DEALLOCATE (eigval)
     442        1260 :       DEALLOCATE (work)
     443        1260 :       DEALLOCATE (dr)
     444             : 
     445             :       CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
     446        1260 :                                         "PRINT%PROGRAM_RUN_INFO")
     447        1260 :       CALL timestop(handle)
     448             : 
     449        5040 :    END SUBROUTINE geoopt_bfgs
     450             : 
     451             : ! **************************************************************************************************
     452             : !> \brief ...
     453             : !> \param ndf ...
     454             : !> \param dg ...
     455             : !> \param eigval ...
     456             : !> \param work ...
     457             : !> \param eigvec_mat ...
     458             : !> \param g ...
     459             : !> \param para_env ...
     460             : ! **************************************************************************************************
     461        3424 :    SUBROUTINE rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
     462             : 
     463             :       INTEGER, INTENT(IN)                                :: ndf
     464             :       REAL(KIND=dp), INTENT(INOUT)                       :: dg(ndf), eigval(ndf), work(ndf)
     465             :       TYPE(cp_fm_type), INTENT(IN)                       :: eigvec_mat
     466             :       REAL(KIND=dp), INTENT(INOUT)                       :: g(ndf)
     467             :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
     468             : 
     469             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'rat_fun_opt'
     470             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp
     471             : 
     472             :       INTEGER                                            :: handle, i, indf, iref, iter, j, k, l, &
     473             :                                                             maxit, ncol_local, nrow_local
     474        1712 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     475             :       LOGICAL                                            :: bisec, fail, set
     476             :       REAL(KIND=dp)                                      :: fun, fun1, fun2, fun3, fung, lam1, lam2, &
     477             :                                                             ln, lp, ssize, step, stol
     478        1712 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER            :: local_data
     479             : 
     480        1712 :       CALL timeset(routineN, handle)
     481             : 
     482        1712 :       stol = 1.0E-8_dp
     483        1712 :       ssize = 0.2_dp
     484        1712 :       maxit = 999
     485        1712 :       fail = .FALSE.
     486        1712 :       bisec = .FALSE.
     487             : 
     488      131300 :       dg = 0._dp
     489             : 
     490             :       CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, &
     491        1712 :                           local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
     492             : 
     493       87716 :       DO i = 1, nrow_local
     494       86004 :          j = row_indices(i)
     495    17038748 :          DO k = 1, ncol_local
     496    16951032 :             l = col_indices(k)
     497    17037036 :             dg(l) = dg(l) + local_data(i, k)*g(j)
     498             :          END DO
     499             :       END DO
     500        1712 :       CALL para_env%sum(dg)
     501             : 
     502        1712 :       set = .FALSE.
     503             : 
     504             :       DO
     505             : 
     506             : !   calculating Lambda
     507             : 
     508        1712 :          lp = 0.0_dp
     509        1712 :          iref = 1
     510        1712 :          ln = 0.0_dp
     511        1712 :          IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp
     512             : 
     513        1712 :          iter = 0
     514             :          DO
     515        5189 :             iter = iter + 1
     516        5189 :             fun = 0.0_dp
     517        5189 :             fung = 0.0_dp
     518      345503 :             DO indf = 1, ndf
     519      340314 :                fun = fun + dg(indf)**2/(ln - eigval(indf))
     520      345503 :                fung = fung - dg(indf)**2/(ln - eigval(indf)**2)
     521             :             END DO
     522        5189 :             fun = fun - ln
     523        5189 :             fung = fung - one
     524        5189 :             step = fun/fung
     525        5189 :             ln = ln - step
     526        5189 :             IF (ABS(step) < stol) GOTO 200
     527        3482 :             IF (iter >= maxit) EXIT
     528             :          END DO
     529             : 100      CONTINUE
     530           5 :          bisec = .TRUE.
     531           5 :          iter = 0
     532           5 :          maxit = 9999
     533           5 :          lam1 = 0.0_dp
     534           5 :          IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp
     535             :          fun1 = 0.0_dp
     536         155 :          DO indf = 1, ndf
     537         155 :             fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf))
     538             :          END DO
     539           5 :          fun1 = fun1 - lam1
     540           5 :          step = ABS(lam1)/1000.0_dp
     541           5 :          IF (step < ssize) step = ssize
     542             :          DO
     543           5 :             iter = iter + 1
     544           5 :             IF (iter > maxit) THEN
     545             :                ln = 0.0_dp
     546        1712 :                lp = 0.0_dp
     547             :                fail = .TRUE.
     548             :                GOTO 300
     549             :             END IF
     550           5 :             fun2 = 0.0_dp
     551           5 :             lam2 = lam1 - iter*step
     552         155 :             DO indf = 1, ndf
     553         155 :                fun2 = fun2 + eigval(indf)**2/(lam2 - eigval(indf))
     554             :             END DO
     555           5 :             fun2 = fun2 - lam2
     556        1717 :             IF (fun2*fun1 < 0.0_dp) THEN
     557             :                iter = 0
     558             :                DO
     559         125 :                   iter = iter + 1
     560         125 :                   IF (iter > maxit) THEN
     561             :                      ln = 0.0_dp
     562             :                      lp = 0.0_dp
     563             :                      fail = .TRUE.
     564             :                      GOTO 300
     565             :                   END IF
     566         125 :                   step = (lam1 + lam2)/2
     567         125 :                   fun3 = 0.0_dp
     568        3875 :                   DO indf = 1, ndf
     569        3875 :                      fun3 = fun3 + dg(indf)**2/(step - eigval(indf))
     570             :                   END DO
     571         125 :                   fun3 = fun3 - step
     572             : 
     573         125 :                   IF (ABS(step - lam2) < stol) THEN
     574             :                      ln = step
     575             :                      GOTO 200
     576             :                   END IF
     577             : 
     578         120 :                   IF (fun3*fun1 < stol) THEN
     579             :                      lam2 = step
     580             :                   ELSE
     581         120 :                      lam1 = step
     582             :                   END IF
     583             :                END DO
     584             :             END IF
     585             :          END DO
     586             : 
     587             : 200      CONTINUE
     588        1717 :          IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. &
     589             :                                        (eigval(iref) > 0.0_dp))) THEN
     590             : 
     591           5 :             IF (.NOT. bisec) GOTO 100
     592             :             ln = 0.0_dp
     593             :             lp = 0.0_dp
     594             :             fail = .TRUE.
     595             :          END IF
     596             : 
     597             : 300      CONTINUE
     598             : 
     599        1712 :          IF (fail .AND. .NOT. set) THEN
     600           0 :             set = .TRUE.
     601           0 :             DO indf = 1, ndf
     602           0 :                eigval(indf) = eigval(indf)*work(indf)
     603             :             END DO
     604             :             CYCLE
     605             :          END IF
     606             : 
     607        1712 :          IF (.NOT. set) THEN
     608      131300 :             work(1:ndf) = one
     609             :          END IF
     610             : 
     611      131300 :          DO indf = 1, ndf
     612      131300 :             eigval(indf) = eigval(indf) - ln
     613             :          END DO
     614             :          EXIT
     615             :       END DO
     616             : 
     617        1712 :       CALL timestop(handle)
     618             : 
     619        1712 :    END SUBROUTINE rat_fun_opt
     620             : 
     621             : ! **************************************************************************************************
     622             : !> \brief ...
     623             : !> \param ndf ...
     624             : !> \param dx ...
     625             : !> \param dg ...
     626             : !> \param hess_mat ...
     627             : !> \param work ...
     628             : !> \param para_env ...
     629             : ! **************************************************************************************************
     630       12214 :    SUBROUTINE bfgs(ndf, dx, dg, hess_mat, work, para_env)
     631             :       INTEGER, INTENT(IN)                                :: ndf
     632             :       REAL(KIND=dp), INTENT(INOUT)                       :: dx(ndf), dg(ndf)
     633             :       TYPE(cp_fm_type), INTENT(IN)                       :: hess_mat
     634             :       REAL(KIND=dp), INTENT(INOUT)                       :: work(ndf)
     635             :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
     636             : 
     637             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'bfgs'
     638             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     639             : 
     640             :       INTEGER                                            :: handle, i, j, k, l, ncol_local, &
     641             :                                                             nrow_local
     642        6107 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     643             :       REAL(KIND=dp)                                      :: DDOT, dxw, gdx
     644        6107 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER            :: local_hes
     645             : 
     646        6107 :       CALL timeset(routineN, handle)
     647             : 
     648             :       CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
     649        6107 :                           local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local)
     650             : 
     651      606722 :       work = zero
     652      341693 :       DO i = 1, nrow_local
     653      335586 :          j = row_indices(i)
     654    54739025 :          DO k = 1, ncol_local
     655    54397332 :             l = col_indices(k)
     656    54732918 :             work(j) = work(j) + local_hes(i, k)*dx(l)
     657             :          END DO
     658             :       END DO
     659             : 
     660        6107 :       CALL para_env%sum(work)
     661             : 
     662        6107 :       gdx = DDOT(ndf, dg, 1, dx, 1)
     663        6107 :       gdx = one/gdx
     664        6107 :       dxw = DDOT(ndf, dx, 1, work, 1)
     665        6107 :       dxw = one/dxw
     666             : 
     667      341693 :       DO i = 1, nrow_local
     668      335586 :          j = row_indices(i)
     669    54739025 :          DO k = 1, ncol_local
     670    54397332 :             l = col_indices(k)
     671             :             local_hes(i, k) = local_hes(i, k) + gdx*dg(j)*dg(l) - &
     672    54732918 :                               dxw*work(j)*work(l)
     673             :          END DO
     674             :       END DO
     675             : 
     676        6107 :       CALL timestop(handle)
     677             : 
     678        6107 :    END SUBROUTINE bfgs
     679             : 
     680             : ! **************************************************************************************************
     681             : !> \brief ...
     682             : !> \param ndf ...
     683             : !> \param eigval ...
     684             : !> \param work ...
     685             : ! **************************************************************************************************
     686        1712 :    SUBROUTINE set_hes_eig(ndf, eigval, work)
     687             :       INTEGER, INTENT(IN)                                :: ndf
     688             :       REAL(KIND=dp), INTENT(INOUT)                       :: eigval(ndf), work(ndf)
     689             : 
     690             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'set_hes_eig'
     691             :       REAL(KIND=dp), PARAMETER                           :: max_neg = -0.5_dp, max_pos = 5.0_dp, &
     692             :                                                             min_eig = 0.005_dp, one = 1.0_dp
     693             : 
     694             :       INTEGER                                            :: handle, indf
     695             :       LOGICAL                                            :: neg
     696             : 
     697        1712 :       CALL timeset(routineN, handle)
     698             : 
     699      131300 :       DO indf = 1, ndf
     700      129588 :          IF (eigval(indf) < 0.0_dp) neg = .TRUE.
     701      131300 :          IF (eigval(indf) > 1000.0_dp) eigval(indf) = 1000.0_dp
     702             :       END DO
     703      131300 :       DO indf = 1, ndf
     704      131300 :          IF (eigval(indf) < 0.0_dp) THEN
     705           5 :             IF (eigval(indf) < max_neg) THEN
     706           1 :                eigval(indf) = max_neg
     707           4 :             ELSE IF (eigval(indf) > -min_eig) THEN
     708           1 :                eigval(indf) = -min_eig
     709             :             END IF
     710      129583 :          ELSE IF (eigval(indf) < 1000.0_dp) THEN
     711      129583 :             IF (eigval(indf) < min_eig) THEN
     712         548 :                eigval(indf) = min_eig
     713      129035 :             ELSE IF (eigval(indf) > max_pos) THEN
     714           0 :                eigval(indf) = max_pos
     715             :             END IF
     716             :          END IF
     717             :       END DO
     718             : 
     719      131300 :       DO indf = 1, ndf
     720      131300 :          IF (eigval(indf) < 0.0_dp) THEN
     721           5 :             work(indf) = -one
     722             :          ELSE
     723      129583 :             work(indf) = one
     724             :          END IF
     725             :       END DO
     726             : 
     727        1712 :       CALL timestop(handle)
     728             : 
     729        1712 :    END SUBROUTINE set_hes_eig
     730             : 
     731             : ! **************************************************************************************************
     732             : !> \brief ...
     733             : !> \param ndf ...
     734             : !> \param eigval ...
     735             : !> \param eigvec_mat ...
     736             : !> \param hess_tmp ...
     737             : !> \param dr ...
     738             : !> \param g ...
     739             : !> \param para_env ...
     740             : !> \param use_rfo ...
     741             : ! **************************************************************************************************
     742       14714 :    SUBROUTINE geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
     743             : 
     744             :       INTEGER, INTENT(IN)                                :: ndf
     745             :       REAL(KIND=dp), INTENT(INOUT)                       :: eigval(ndf)
     746             :       TYPE(cp_fm_type), INTENT(IN)                       :: eigvec_mat, hess_tmp
     747             :       REAL(KIND=dp), INTENT(INOUT)                       :: dr(ndf), g(ndf)
     748             :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
     749             :       LOGICAL                                            :: use_rfo
     750             : 
     751             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     752             : 
     753             :       INTEGER                                            :: i, indf, j, k, l, ncol_local, nrow_local
     754        7357 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     755        7357 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER            :: local_data
     756             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     757             :       TYPE(cp_fm_type)                          :: tmp
     758             : 
     759        7357 :       CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
     760        7357 :       IF (use_rfo) THEN
     761      131300 :          DO indf = 1, ndf
     762      131300 :             eigval(indf) = one/eigval(indf)
     763             :          END DO
     764             :       ELSE
     765      564968 :          DO indf = 1, ndf
     766      564968 :             eigval(indf) = one/MAX(0.0001_dp, eigval(indf))
     767             :          END DO
     768             :       END IF
     769             : 
     770        7357 :       CALL cp_fm_column_scale(hess_tmp, eigval)
     771        7357 :       CALL cp_fm_get_info(eigvec_mat, matrix_struct=matrix_struct)
     772        7357 :       CALL cp_fm_create(tmp, matrix_struct, name="tmp")
     773             : 
     774        7357 :       CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp)
     775             : 
     776        7357 :       CALL cp_fm_transpose(tmp, hess_tmp)
     777        7357 :       CALL cp_fm_release(tmp)
     778             : 
     779             : !    ** New step **
     780             : 
     781             :       CALL cp_fm_get_info(hess_tmp, row_indices=row_indices, col_indices=col_indices, &
     782        7357 :                           local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
     783             : 
     784      696268 :       dr = 0.0_dp
     785      389521 :       DO i = 1, nrow_local
     786      382164 :          j = row_indices(i)
     787    64565677 :          DO k = 1, ncol_local
     788    64176156 :             l = col_indices(k)
     789    64558320 :             dr(j) = dr(j) - local_data(i, k)*g(l)
     790             :          END DO
     791             :       END DO
     792             : 
     793        7357 :       CALL para_env%sum(dr)
     794             : 
     795        7357 :    END SUBROUTINE geoopt_get_step
     796             : 
     797             : ! **************************************************************************************************
     798             : !> \brief ...
     799             : !> \param ndf ...
     800             : !> \param step ...
     801             : !> \param rad ...
     802             : !> \param rat ...
     803             : !> \param dr ...
     804             : !> \param output_unit ...
     805             : ! **************************************************************************************************
     806        7357 :    SUBROUTINE trust_radius(ndf, step, rad, rat, dr, output_unit)
     807             :       INTEGER, INTENT(IN)                                :: ndf
     808             :       REAL(KIND=dp), INTENT(INOUT)                       :: step, rad, rat, dr(ndf)
     809             :       INTEGER, INTENT(IN)                                :: output_unit
     810             : 
     811             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'trust_radius'
     812             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp
     813             : 
     814             :       INTEGER                                            :: handle
     815             :       REAL(KIND=dp)                                      :: scal
     816             : 
     817        7357 :       CALL timeset(routineN, handle)
     818             : 
     819      703625 :       step = MAXVAL(ABS(dr))
     820        7357 :       scal = MAX(one, rad/step)
     821             : 
     822        7357 :       IF (step > rad) THEN
     823         465 :          rat = rad/step
     824         465 :          CALL DSCAL(ndf, rat, dr, 1)
     825         465 :          step = rad
     826         465 :          IF (output_unit > 0) THEN
     827             :             WRITE (unit=output_unit, FMT="(/,T2,A,F8.5)") &
     828         219 :                " Step is scaled; Scaling factor = ", rat
     829         219 :             CALL m_flush(output_unit)
     830             :          END IF
     831             :       END IF
     832        7357 :       CALL timestop(handle)
     833             : 
     834        7357 :    END SUBROUTINE trust_radius
     835             : 
     836             : ! **************************************************************************************************
     837             : !> \brief ...
     838             : !> \param ndf ...
     839             : !> \param work ...
     840             : !> \param hess_mat ...
     841             : !> \param dr ...
     842             : !> \param g ...
     843             : !> \param conv ...
     844             : !> \param pred ...
     845             : !> \param para_env ...
     846             : ! **************************************************************************************************
     847       14714 :    SUBROUTINE energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
     848             : 
     849             :       INTEGER, INTENT(IN)                                :: ndf
     850             :       REAL(KIND=dp), INTENT(INOUT)                       :: work(ndf)
     851             :       TYPE(cp_fm_type), INTENT(IN)                           :: hess_mat
     852             :       REAL(KIND=dp), INTENT(INOUT)                       :: dr(ndf), g(ndf)
     853             :       LOGICAL, INTENT(INOUT)                             :: conv
     854             :       REAL(KIND=dp), INTENT(INOUT)                       :: pred
     855             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     856             : 
     857             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'energy_predict'
     858             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
     859             : 
     860             :       INTEGER                                            :: handle, i, j, k, l, ncol_local, &
     861             :                                                             nrow_local
     862        7357 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     863             :       REAL(KIND=dp)                                      :: DDOT, ener1, ener2
     864        7357 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER            :: local_data
     865             : 
     866        7357 :       CALL timeset(routineN, handle)
     867             : 
     868        7357 :       ener1 = DDOT(ndf, g, 1, dr, 1)
     869             : 
     870             :       CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
     871        7357 :                           local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
     872             : 
     873      696268 :       work = zero
     874      389521 :       DO i = 1, nrow_local
     875      382164 :          j = row_indices(i)
     876    64565677 :          DO k = 1, ncol_local
     877    64176156 :             l = col_indices(k)
     878    64558320 :             work(j) = work(j) + local_data(i, k)*dr(l)
     879             :          END DO
     880             :       END DO
     881             : 
     882        7357 :       CALL para_env%sum(work)
     883        7357 :       ener2 = DDOT(ndf, dr, 1, work, 1)
     884        7357 :       pred = ener1 + 0.5_dp*ener2
     885        7357 :       conv = .FALSE.
     886        7357 :       CALL timestop(handle)
     887             : 
     888        7357 :    END SUBROUTINE energy_predict
     889             : 
     890             : ! **************************************************************************************************
     891             : !> \brief ...
     892             : !> \param rat ...
     893             : !> \param rad ...
     894             : !> \param step ...
     895             : !> \param ediff ...
     896             : ! **************************************************************************************************
     897        1582 :    SUBROUTINE update_trust_rad(rat, rad, step, ediff)
     898             : 
     899             :       REAL(KIND=dp), INTENT(INOUT)                       :: rat, rad, step, ediff
     900             : 
     901             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'update_trust_rad'
     902             :       REAL(KIND=dp), PARAMETER                           :: max_trust = 1.0_dp, min_trust = 0.1_dp
     903             : 
     904             :       INTEGER                                            :: handle
     905             : 
     906        1582 :       CALL timeset(routineN, handle)
     907             : 
     908        1582 :       IF (rat > 4.0_dp) THEN
     909           0 :          IF (ediff < 0.0_dp) THEN
     910           0 :             rad = step*0.5_dp
     911             :          ELSE
     912           0 :             rad = step*0.25_dp
     913             :          END IF
     914        1582 :       ELSE IF (rat > 2.0_dp) THEN
     915           0 :          IF (ediff < 0.0_dp) THEN
     916           0 :             rad = step*0.75_dp
     917             :          ELSE
     918           0 :             rad = step*0.5_dp
     919             :          END IF
     920        1582 :       ELSE IF (rat > 4.0_dp/3.0_dp) THEN
     921           0 :          IF (ediff < 0.0_dp) THEN
     922           0 :             rad = step
     923             :          ELSE
     924           0 :             rad = step*0.75_dp
     925             :          END IF
     926        1582 :       ELSE IF (rat > 10.0_dp/9.0_dp) THEN
     927           0 :          IF (ediff < 0.0_dp) THEN
     928           0 :             rad = step*1.25_dp
     929             :          ELSE
     930           0 :             rad = step
     931             :          END IF
     932        1582 :       ELSE IF (rat > 0.9_dp) THEN
     933          98 :          IF (ediff < 0.0_dp) THEN
     934          97 :             rad = step*1.5_dp
     935             :          ELSE
     936           1 :             rad = step*1.25_dp
     937             :          END IF
     938        1484 :       ELSE IF (rat > 0.75_dp) THEN
     939         200 :          IF (ediff < 0.0_dp) THEN
     940         197 :             rad = step*1.25_dp
     941             :          ELSE
     942           3 :             rad = step
     943             :          END IF
     944        1284 :       ELSE IF (rat > 0.5_dp) THEN
     945         285 :          IF (ediff < 0.0_dp) THEN
     946         280 :             rad = step
     947             :          ELSE
     948           5 :             rad = step*0.75_dp
     949             :          END IF
     950         999 :       ELSE IF (rat > 0.25_dp) THEN
     951           7 :          IF (ediff < 0.0_dp) THEN
     952           7 :             rad = step*0.75_dp
     953             :          ELSE
     954           0 :             rad = step*0.5_dp
     955             :          END IF
     956         992 :       ELSE IF (ediff < 0.0_dp) THEN
     957         992 :          rad = step*0.5_dp
     958             :       ELSE
     959           0 :          rad = step*0.25_dp
     960             :       END IF
     961             : 
     962        1582 :       rad = MAX(rad, min_trust)
     963        1582 :       rad = MIN(rad, max_trust)
     964        1582 :       CALL timestop(handle)
     965             : 
     966        1582 :    END SUBROUTINE update_trust_rad
     967             : 
     968             : ! **************************************************************************************************
     969             : 
     970             : ! **************************************************************************************************
     971             : !> \brief ...
     972             : !> \param geo_section ...
     973             : !> \param hess_mat ...
     974             : !> \param logger ...
     975             : ! **************************************************************************************************
     976        5023 :    SUBROUTINE write_bfgs_hessian(geo_section, hess_mat, logger)
     977             : 
     978             :       TYPE(section_vals_type), POINTER                   :: geo_section
     979             :       TYPE(cp_fm_type), INTENT(IN)                          :: hess_mat
     980             :       TYPE(cp_logger_type), POINTER                      :: logger
     981             : 
     982             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bfgs_hessian'
     983             : 
     984             :       INTEGER                                            :: handle, hesunit
     985             : 
     986        5023 :       CALL timeset(routineN, handle)
     987             : 
     988             :       hesunit = cp_print_key_unit_nr(logger, geo_section, "BFGS%RESTART", &
     989             :                                      extension=".Hessian", file_form="UNFORMATTED", file_action="WRITE", &
     990        5023 :                                      file_position="REWIND")
     991             : 
     992        5023 :       CALL cp_fm_write_unformatted(hess_mat, hesunit)
     993             : 
     994        5023 :       CALL cp_print_key_finished_output(hesunit, logger, geo_section, "BFGS%RESTART")
     995             : 
     996        5023 :       CALL timestop(handle)
     997             : 
     998        5023 :    END SUBROUTINE write_bfgs_hessian
     999             : ! **************************************************************************************************
    1000             : 
    1001             : ! **************************************************************************************************
    1002             : !> \brief ...
    1003             : !> \param force_env ...
    1004             : !> \param hess_mat ...
    1005             : ! **************************************************************************************************
    1006        1048 :    SUBROUTINE construct_initial_hess(force_env, hess_mat)
    1007             : 
    1008             :       TYPE(force_env_type), POINTER                      :: force_env
    1009             :       TYPE(cp_fm_type), INTENT(IN)                          :: hess_mat
    1010             : 
    1011             :       INTEGER                                            :: i, iat_col, iat_row, iglobal, iind, j, &
    1012             :                                                             jat_row, jglobal, jind, k, natom, &
    1013             :                                                             ncol_local, nrow_local, z
    1014        1048 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: at_row
    1015        1048 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1016        1048 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: d_ij, rho_ij
    1017        1048 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: r_ij
    1018             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: alpha, r0
    1019        1048 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER            :: fixed, local_data
    1020             :       TYPE(cell_type), POINTER                           :: cell
    1021             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1022             :       TYPE(particle_list_type), POINTER                  :: particles
    1023             : 
    1024        1048 :       CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
    1025             :       CALL cp_subsys_get(subsys, &
    1026        1048 :                          particles=particles)
    1027             : 
    1028        4192 :       alpha(1, :) = (/1._dp, 0.3949_dp, 0.3949_dp/)
    1029        4192 :       alpha(2, :) = (/0.3494_dp, 0.2800_dp, 0.2800_dp/)
    1030        4192 :       alpha(3, :) = (/0.3494_dp, 0.2800_dp, 0.1800_dp/)
    1031             : 
    1032        4192 :       r0(1, :) = (/1.35_dp, 2.10_dp, 2.53_dp/)
    1033        4192 :       r0(2, :) = (/2.10_dp, 2.87_dp, 3.40_dp/)
    1034        4192 :       r0(3, :) = (/2.53_dp, 3.40_dp, 3.40_dp/)
    1035             : 
    1036             :       CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
    1037        1048 :                           local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
    1038        1048 :       natom = particles%n_els
    1039        3144 :       ALLOCATE (at_row(natom))
    1040        4192 :       ALLOCATE (rho_ij(natom, natom))
    1041        4192 :       ALLOCATE (d_ij(natom, natom))
    1042        5240 :       ALLOCATE (r_ij(natom, natom, 3))
    1043        3144 :       ALLOCATE (fixed(3, natom))
    1044       74080 :       fixed = 1.0_dp
    1045        1048 :       CALL fix_atom_control(force_env, fixed)
    1046        4192 :       DO i = 1, 3
    1047      113740 :          CALL hess_mat%matrix_struct%para_env%min(fixed(i, :))
    1048             :       END DO
    1049      920272 :       rho_ij = 0
    1050             :       !XXXX insert proper rows !XXX
    1051       19306 :       at_row = 3
    1052       19306 :       DO i = 1, natom
    1053       18258 :          CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, z=z)
    1054       18258 :          IF (z .LE. 10) at_row(i) = 2
    1055       37564 :          IF (z .LE. 2) at_row(i) = 1
    1056             :       END DO
    1057       18258 :       DO i = 2, natom
    1058       17210 :          iat_row = at_row(i)
    1059      459612 :          DO j = 1, i - 1
    1060      441354 :             jat_row = at_row(j)
    1061             :             !pbc for a distance vector
    1062     1765416 :             r_ij(j, i, :) = pbc(particles%els(i)%r, particles%els(j)%r, cell)
    1063     1765416 :             r_ij(i, j, :) = -r_ij(j, i, :)
    1064     1765416 :             d_ij(j, i) = SQRT(DOT_PRODUCT(r_ij(j, i, :), r_ij(j, i, :)))
    1065      441354 :             d_ij(i, j) = d_ij(j, i)
    1066      441354 :             rho_ij(j, i) = EXP(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2))
    1067      458564 :             rho_ij(i, j) = rho_ij(j, i)
    1068             :          END DO
    1069             :       END DO
    1070       55822 :       DO i = 1, ncol_local
    1071       54774 :          iglobal = col_indices(i)
    1072       54774 :          iind = MOD(iglobal - 1, 3) + 1
    1073       54774 :          iat_col = (iglobal + 2)/3
    1074       54774 :          IF (iat_col .GT. natom) CYCLE
    1075     4215199 :          DO j = 1, nrow_local
    1076     4159377 :             jglobal = row_indices(j)
    1077     4159377 :             jind = MOD(jglobal - 1, 3) + 1
    1078     4159377 :             iat_row = (jglobal + 2)/3
    1079     4159377 :             IF (iat_row .GT. natom) CYCLE
    1080     4159377 :             IF (iat_row .NE. iat_col) THEN
    1081     4069926 :                IF (d_ij(iat_row, iat_col) .LT. 6.0_dp) &
    1082             :                   local_data(j, i) = local_data(j, i) + &
    1083      372492 :                                      angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
    1084             :             ELSE
    1085             :                local_data(j, i) = local_data(j, i) + &
    1086       89451 :                                   angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
    1087             :             END IF
    1088     4159377 :             IF (iat_col .NE. iat_row) THEN
    1089     4069926 :                IF (d_ij(iat_row, iat_col) .LT. 6.0_dp) &
    1090             :                   local_data(j, i) = local_data(j, i) - &
    1091             :                                      dist_second_deriv(r_ij(iat_col, iat_row, :), &
    1092     2607444 :                                                        iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col))
    1093             :             ELSE
    1094     4248828 :                DO k = 1, natom
    1095     4159377 :                   IF (k == iat_col) CYCLE
    1096     4069926 :                   IF (d_ij(iat_row, k) .LT. 6.0_dp) &
    1097             :                      local_data(j, i) = local_data(j, i) + &
    1098             :                                         dist_second_deriv(r_ij(iat_col, k, :), &
    1099     2696895 :                                                           iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k))
    1100             :                END DO
    1101             :             END IF
    1102     4214151 :             IF (fixed(jind, iat_row) .LT. 0.5_dp .OR. fixed(iind, iat_col) .LT. 0.5_dp) THEN
    1103       10161 :                local_data(j, i) = 0.0_dp
    1104       10161 :                IF (jind == iind .AND. iat_row == iat_col) local_data(j, i) = 1.0_dp
    1105             :             END IF
    1106             :          END DO
    1107             :       END DO
    1108        1048 :       DEALLOCATE (fixed)
    1109        1048 :       DEALLOCATE (rho_ij)
    1110        1048 :       DEALLOCATE (d_ij)
    1111        1048 :       DEALLOCATE (r_ij)
    1112        1048 :       DEALLOCATE (at_row)
    1113             : 
    1114        3144 :    END SUBROUTINE construct_initial_hess
    1115             : 
    1116             : ! **************************************************************************************************
    1117             : !> \brief ...
    1118             : !> \param r1 ...
    1119             : !> \param i ...
    1120             : !> \param j ...
    1121             : !> \param d ...
    1122             : !> \param rho ...
    1123             : !> \return ...
    1124             : ! **************************************************************************************************
    1125      744984 :    FUNCTION dist_second_deriv(r1, i, j, d, rho) RESULT(deriv)
    1126             :       REAL(KIND=dp), DIMENSION(3)                        :: r1
    1127             :       INTEGER                                            :: i, j
    1128             :       REAL(KIND=dp)                                      :: d, rho, deriv
    1129             : 
    1130      744984 :       deriv = 0.45_dp*rho*(r1(i)*r1(j))/d**2
    1131      744984 :    END FUNCTION
    1132             : 
    1133             : ! **************************************************************************************************
    1134             : !> \brief ...
    1135             : !> \param r_ij ...
    1136             : !> \param d_ij ...
    1137             : !> \param rho_ij ...
    1138             : !> \param idir ...
    1139             : !> \param jdir ...
    1140             : !> \param iat_der ...
    1141             : !> \param jat_der ...
    1142             : !> \param natom ...
    1143             : !> \return ...
    1144             : ! **************************************************************************************************
    1145      461943 :    FUNCTION angle_second_deriv(r_ij, d_ij, rho_ij, idir, jdir, iat_der, jat_der, natom) RESULT(deriv)
    1146             :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: r_ij
    1147             :       REAL(KIND=dp), DIMENSION(:, :)                     :: d_ij, rho_ij
    1148             :       INTEGER                                            :: idir, jdir, iat_der, jat_der, natom
    1149             :       REAL(KIND=dp)                                      :: deriv
    1150             : 
    1151             :       INTEGER                                            :: i, iat, idr, j, jat, jdr
    1152             :       REAL(KIND=dp)                                      :: d12, d23, d31, D_mat(3, 2), denom1, &
    1153             :                                                             denom2, denom3, ka1, ka2, ka3, rho12, &
    1154             :                                                             rho23, rho31, rsst1, rsst2, rsst3
    1155             :       REAL(KIND=dp), DIMENSION(3)                        :: r12, r23, r31
    1156             : 
    1157      461943 :       deriv = 0._dp
    1158      461943 :       IF (iat_der == jat_der) THEN
    1159     4159377 :          DO i = 1, natom - 1
    1160     4069926 :             IF (rho_ij(iat_der, i) .LT. 0.00001) CYCLE
    1161    25833726 :             DO j = i + 1, natom
    1162    24784281 :                IF (rho_ij(iat_der, j) .LT. 0.00001) CYCLE
    1163     6121404 :                IF (i == iat_der .OR. j == iat_der) CYCLE
    1164     6121404 :                IF (iat_der .LT. i .OR. iat_der .GT. j) THEN
    1165    39396420 :                   r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :)
    1166     3939642 :                   d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der)
    1167     3939642 :                   rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der)
    1168             :                ELSE
    1169    21817620 :                   r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :)
    1170     2181762 :                   d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der)
    1171     2181762 :                   rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der)
    1172             :                END IF
    1173     6121404 :                ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
    1174    61214040 :                rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
    1175     6121404 :                denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
    1176     6121404 :                denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
    1177     6121404 :                denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
    1178     6121404 :                denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
    1179     6121404 :                denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
    1180     6121404 :                D_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23)
    1181     6121404 :                D_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23)
    1182     6121404 :                D_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3)
    1183     6121404 :                D_mat(2, 2) = -r23(jdir)/(d23*d31) + rsst2*r31(jdir)/(d23*d31**3)
    1184             :                D_mat(3, 1) = (r31(idir) - r12(idir))/(d31*d12) + rsst3*r31(idir)/(d31**3*d12) - &
    1185     6121404 :                              rsst3*r12(idir)/(d31*d12**3)
    1186             :                D_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - &
    1187     6121404 :                              rsst3*r12(jdir)/(d31*d12**3)
    1188     6121404 :                IF (ABS(denom1) .LE. 0.011_dp) D_mat(1, 1) = 0.0_dp
    1189     6121404 :                IF (ABS(denom2) .LE. 0.011_dp) D_mat(2, 1) = 0.0_dp
    1190     6121404 :                IF (ABS(denom3) .LE. 0.011_dp) D_mat(3, 1) = 0.0_dp
    1191             :                deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
    1192             :                        ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
    1193    28854207 :                        ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
    1194             : 
    1195             :             END DO
    1196             :          END DO
    1197             :       ELSE
    1198    10621998 :          DO i = 1, natom
    1199    10249506 :             IF (i == iat_der .OR. i == jat_der) CYCLE
    1200     9504522 :             IF (jat_der .LT. iat_der) THEN
    1201     4752261 :                iat = jat_der; jat = iat_der; idr = jdir; jdr = idir
    1202             :             ELSE
    1203     4752261 :                iat = iat_der; jat = jat_der; idr = idir; jdr = jdir
    1204             :             END IF
    1205     9504522 :             IF (jat .LT. i .OR. iat .GT. i) THEN
    1206    71616240 :                r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :)
    1207     7161624 :                d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat)
    1208     7161624 :                rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat)
    1209             :             ELSE
    1210    23428980 :                r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :)
    1211     2342898 :                d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat)
    1212     2342898 :                rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat)
    1213             :             END IF
    1214     9504522 :             ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
    1215    95045220 :             rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
    1216     9504522 :             denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
    1217     9504522 :             denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
    1218     9504522 :             denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
    1219     9504522 :             denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
    1220     9504522 :             denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
    1221     9504522 :             D_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23)
    1222     9504522 :             D_mat(2, 1) = -r23(idr)/(d23*d31) + rsst2*r31(idr)/(d23*d31**3)
    1223             :             D_mat(3, 1) = (r31(idr) - r12(idr))/(d31*d12) + rsst3*r31(idr)/(d31**3*d12) - &
    1224     9504522 :                           rsst3*r12(idr)/(d31*d12**3)
    1225     9504522 :             IF (jat .LT. i .OR. iat .GT. i) THEN
    1226             :                D_mat(1, 2) = (r12(jdr) - r23(jdr))/(d12*d23) + rsst1*r12(jdr)/(d12**3*d23) - &
    1227     7161624 :                              rsst1*r23(jdr)/(d12*d23**3)
    1228     7161624 :                D_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31)
    1229     7161624 :                D_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3)
    1230             :             ELSE
    1231     2342898 :                D_mat(1, 2) = -r12(jdr)/(d12*d23) + rsst1*r23(jdr)/(d12*d23**3)
    1232             :                D_mat(2, 2) = (r23(jdr) - r31(jdr))/(d23*d31) + rsst2*r23(jdr)/(d23**3*d31) - &
    1233     2342898 :                              rsst2*r31(jdr)/(d23*d31**3)
    1234     2342898 :                D_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12)
    1235             :             END IF
    1236     9504522 :             IF (ABS(denom1) .LE. 0.011_dp) D_mat(1, 1) = 0.0_dp
    1237     9504522 :             IF (ABS(denom2) .LE. 0.011_dp) D_mat(2, 1) = 0.0_dp
    1238     9504522 :             IF (ABS(denom3) .LE. 0.011_dp) D_mat(3, 1) = 0.0_dp
    1239             : 
    1240             :             deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
    1241             :                     ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
    1242    10621998 :                     ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
    1243             :          END DO
    1244             :       END IF
    1245      461943 :       deriv = 0.25_dp*deriv
    1246             : 
    1247      461943 :    END FUNCTION angle_second_deriv
    1248             : 
    1249             : END MODULE bfgs_optimizer

Generated by: LCOV version 1.15