LCOV - code coverage report
Current view: top level - src/motion - cg_optimizer.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 107 116 92.2 %
Date: 2024-04-20 06:29:22 Functions: 2 2 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  Conjugate Gradients
      10             : !> \author Teodoro Laino [teo]
      11             : !>      10.2005
      12             : ! **************************************************************************************************
      13             : MODULE cg_optimizer
      14             : 
      15             :    USE cell_types, ONLY: cell_type
      16             :    USE cg_utils, ONLY: cg_linmin, &
      17             :                        get_conjugate_direction
      18             :    USE cp_external_control, ONLY: external_control
      19             :    USE cp_log_handling, ONLY: cp_get_default_logger, &
      20             :                               cp_logger_type
      21             :    USE cp_output_handling, ONLY: cp_iterate, &
      22             :                                  cp_print_key_finished_output, &
      23             :                                  cp_print_key_unit_nr
      24             :    USE cp_subsys_types, ONLY: cp_subsys_type
      25             :    USE force_env_types, ONLY: force_env_get, &
      26             :                               force_env_type
      27             :    USE global_types, ONLY: global_environment_type
      28             :    USE gopt_f_methods, ONLY: gopt_f_ii, &
      29             :                              gopt_f_io, &
      30             :                              gopt_f_io_finalize, &
      31             :                              gopt_f_io_init, &
      32             :                              print_geo_opt_header, &
      33             :                              print_geo_opt_nc
      34             :    USE gopt_f_types, ONLY: gopt_f_type
      35             :    USE gopt_param_types, ONLY: gopt_param_type
      36             :    USE input_constants, ONLY: default_cell_direct_id, &
      37             :                               default_cell_geo_opt_id, &
      38             :                               default_cell_md_id, &
      39             :                               default_cell_method_id, &
      40             :                               default_minimization_method_id, &
      41             :                               default_ts_method_id
      42             :    USE input_section_types, ONLY: section_vals_type, &
      43             :                                   section_vals_val_get, &
      44             :                                   section_vals_val_set
      45             :    USE kinds, ONLY: dp
      46             :    USE machine, ONLY: m_walltime
      47             :    USE message_passing, ONLY: mp_para_env_type
      48             :    USE space_groups, ONLY: identify_space_group, &
      49             :                            print_spgr, &
      50             :                            spgr_apply_rotations_coord, &
      51             :                            spgr_apply_rotations_force
      52             :    USE space_groups_types, ONLY: spgr_type
      53             : #include "../base/base_uses.f90"
      54             : 
      55             :    IMPLICIT NONE
      56             :    PRIVATE
      57             : 
      58             :    #:include "gopt_f77_methods.fypp"
      59             : 
      60             :    PUBLIC :: geoopt_cg
      61             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      62             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cg_optimizer'
      63             : 
      64             : CONTAINS
      65             : 
      66             : ! **************************************************************************************************
      67             : !> \brief Driver for conjugate gradient optimization technique
      68             : !> \param force_env ...
      69             : !> \param gopt_param ...
      70             : !> \param globenv ...
      71             : !> \param geo_section ...
      72             : !> \param gopt_env ...
      73             : !> \param x0 ...
      74             : !> \param do_update ...
      75             : !> \par History
      76             : !>      10.2005 created [tlaino]
      77             : !> \author Teodoro Laino
      78             : ! **************************************************************************************************
      79         956 :    RECURSIVE SUBROUTINE geoopt_cg(force_env, gopt_param, globenv, geo_section, &
      80             :                                   gopt_env, x0, do_update)
      81             : 
      82             :       TYPE(force_env_type), POINTER                      :: force_env
      83             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
      84             :       TYPE(global_environment_type), POINTER             :: globenv
      85             :       TYPE(section_vals_type), POINTER                   :: geo_section
      86             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
      87             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
      88             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: do_update
      89             : 
      90             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'geoopt_cg'
      91             : 
      92             :       INTEGER                                            :: handle, output_unit
      93             :       LOGICAL                                            :: my_do_update
      94             :       TYPE(cp_logger_type), POINTER                      :: logger
      95             :       TYPE(cp_subsys_type), POINTER                      :: subsys
      96             :       TYPE(spgr_type), POINTER                           :: spgr
      97             : 
      98         478 :       CALL timeset(routineN, handle)
      99             : 
     100         478 :       NULLIFY (spgr)
     101         478 :       logger => cp_get_default_logger()
     102         478 :       spgr => gopt_env%spgr
     103             : 
     104             :       output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", &
     105         478 :                                          extension=".geoLog")
     106         478 :       CALL print_geo_opt_header(gopt_env, output_unit, "CONJUGATE GRADIENTS")
     107             : 
     108             :       ! find space_group
     109         478 :       CALL force_env_get(force_env, subsys=subsys)
     110         478 :       CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
     111         478 :       IF (spgr%keep_space_group) THEN
     112           2 :          SELECT CASE (gopt_env%type_id)
     113             :          CASE (default_minimization_method_id, default_ts_method_id)
     114           0 :             CALL force_env_get(force_env, subsys=subsys)
     115           0 :             CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
     116           0 :             CALL spgr_apply_rotations_coord(spgr, x0)
     117           0 :             CALL print_spgr(spgr)
     118             :          CASE (default_cell_method_id)
     119           4 :             SELECT CASE (gopt_env%cell_method_id)
     120             :             CASE (default_cell_direct_id)
     121           2 :                CALL force_env_get(force_env, subsys=subsys)
     122           2 :                CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
     123           2 :                CALL spgr_apply_rotations_coord(spgr, x0)
     124           2 :                CALL print_spgr(spgr)
     125             :             CASE (default_cell_geo_opt_id)
     126           0 :                spgr%keep_space_group = .FALSE.
     127             :             CASE (default_cell_md_id)
     128           0 :                CPABORT("KEEP_SPACE_GROUP not implemented for motion method MD.")
     129             :             CASE DEFAULT
     130           2 :                spgr%keep_space_group = .FALSE.
     131             :             END SELECT
     132             :          CASE DEFAULT
     133           2 :             spgr%keep_space_group = .FALSE.
     134             :          END SELECT
     135             :       END IF
     136             : 
     137             :       CALL cp_cg_main(force_env, x0, gopt_param, output_unit, globenv, &
     138         478 :                       gopt_env, do_update=my_do_update)
     139             : 
     140             :       CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
     141         478 :                                         "PRINT%PROGRAM_RUN_INFO")
     142         478 :       IF (PRESENT(do_update)) do_update = my_do_update
     143             : 
     144         478 :       CALL timestop(handle)
     145             : 
     146         478 :    END SUBROUTINE geoopt_cg
     147             : 
     148             : ! **************************************************************************************************
     149             : !> \brief This really performs the conjugate gradients optimization
     150             : !> \param force_env ...
     151             : !> \param x0 ...
     152             : !> \param gopt_param ...
     153             : !> \param output_unit ...
     154             : !> \param globenv ...
     155             : !> \param gopt_env ...
     156             : !> \param do_update ...
     157             : !> \par History
     158             : !>      10.2005 created [tlaino]
     159             : !> \author Teodoro Laino
     160             : ! **************************************************************************************************
     161         478 :    RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv, &
     162             :                                    gopt_env, do_update)
     163             :       TYPE(force_env_type), POINTER                      :: force_env
     164             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
     165             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     166             :       INTEGER, INTENT(IN)                                :: output_unit
     167             :       TYPE(global_environment_type), POINTER             :: globenv
     168             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     169             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: do_update
     170             : 
     171             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_cg_main'
     172             : 
     173             :       CHARACTER(LEN=5)                                   :: wildcard
     174             :       INTEGER                                            :: handle, iter_nr, its, max_steep_steps, &
     175             :                                                             maxiter
     176             :       LOGICAL                                            :: conv, Fletcher_Reeves, &
     177             :                                                             save_consistent_energy_force, &
     178             :                                                             should_stop
     179             :       REAL(KIND=dp)                                      :: emin, eold, opt_energy, res_lim, t_diff, &
     180             :                                                             t_now, t_old
     181         478 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: xold
     182         478 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: g, h, xi
     183             :       TYPE(cell_type), POINTER                           :: cell
     184             :       TYPE(cp_logger_type), POINTER                      :: logger
     185             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     186             :       TYPE(section_vals_type), POINTER                   :: root_section
     187             :       TYPE(spgr_type), POINTER                           :: spgr
     188             : 
     189         478 :       CALL timeset(routineN, handle)
     190         478 :       t_old = m_walltime()
     191         478 :       NULLIFY (logger, g, h, xi, spgr)
     192         478 :       root_section => force_env%root_section
     193         478 :       logger => cp_get_default_logger()
     194         478 :       conv = .FALSE.
     195         478 :       maxiter = gopt_param%max_iter
     196         478 :       max_steep_steps = gopt_param%max_steep_steps
     197         478 :       Fletcher_Reeves = gopt_param%Fletcher_Reeves
     198         478 :       res_lim = gopt_param%restart_limit
     199        1434 :       ALLOCATE (g(SIZE(x0)))
     200        1434 :       ALLOCATE (h(SIZE(x0)))
     201        1434 :       ALLOCATE (xi(SIZE(x0)))
     202        1434 :       ALLOCATE (xold(SIZE(x0)))
     203         478 :       CALL force_env_get(force_env, cell=cell, subsys=subsys)
     204             : 
     205         478 :       spgr => gopt_env%spgr
     206             :       ! applies rotation matrices to coordinates
     207         478 :       IF (spgr%keep_space_group) THEN
     208           2 :          CALL spgr_apply_rotations_coord(spgr, x0)
     209             :       END IF
     210             : 
     211             :       ! Evaluate energy and forces at the first step
     212             :       ![NB] consistent energies and forces not required for CG, but some line minimizers might set it
     213         478 :       save_consistent_energy_force = gopt_env%require_consistent_energy_force
     214         478 :       gopt_env%require_consistent_energy_force = .FALSE.
     215             : 
     216             :       CALL cp_eval_at(gopt_env, x0, opt_energy, xi, gopt_env%force_env%para_env%mepos, &
     217         478 :                       .FALSE., gopt_env%force_env%para_env)
     218             : 
     219         478 :       gopt_env%require_consistent_energy_force = save_consistent_energy_force
     220             : 
     221             :       ! Symmetrize coordinates and forces
     222         478 :       IF (spgr%keep_space_group) THEN
     223           2 :          CALL spgr_apply_rotations_coord(spgr, x0)
     224           2 :          CALL spgr_apply_rotations_force(spgr, xi)
     225             :       END IF
     226             : 
     227      185012 :       g = -xi
     228      185012 :       h = g
     229      185012 :       xi = h
     230         478 :       emin = HUGE(0.0_dp)
     231         478 :       CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
     232             :       ! Main Loop
     233         478 :       wildcard = "   SD"
     234         478 :       t_now = m_walltime()
     235         478 :       t_diff = t_now - t_old
     236         478 :       t_old = t_now
     237         478 :       CALL gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, used_time=t_diff, its=iter_nr)
     238         478 :       eold = opt_energy
     239        2526 :       DO its = iter_nr + 1, maxiter
     240        2526 :          CALL cp_iterate(logger%iter_info, last=(its == maxiter))
     241        2526 :          CALL section_vals_val_set(gopt_env%geo_section, "STEP_START_VAL", i_val=its)
     242        2526 :          CALL gopt_f_ii(its, output_unit)
     243             : 
     244             :          ! Symmetrize coordinates and forces
     245        2526 :          IF (spgr%keep_space_group) THEN
     246           2 :             CALL spgr_apply_rotations_coord(spgr, x0)
     247           2 :             CALL spgr_apply_rotations_force(spgr, g)
     248           2 :             CALL spgr_apply_rotations_force(spgr, xi)
     249             :          END IF
     250             : 
     251      415038 :          xold(:) = x0
     252             : 
     253             :          ! Line minimization
     254        2526 :          CALL cg_linmin(gopt_env, x0, xi, g, opt_energy, output_unit, gopt_param, globenv)
     255             : 
     256             :          ! Applies rotation matrices to coordinates
     257        2526 :          IF (spgr%keep_space_group) THEN
     258           2 :             CALL spgr_apply_rotations_coord(spgr, x0)
     259             :          END IF
     260             : 
     261             :          ! Check for an external exit command
     262        2526 :          CALL external_control(should_stop, "GEO", globenv=globenv)
     263        2526 :          IF (should_stop) EXIT
     264             : 
     265             :          ! Some IO and Convergence check
     266        2526 :          t_now = m_walltime()
     267        2526 :          t_diff = t_now - t_old
     268        2526 :          t_old = t_now
     269             :          CALL gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
     270             :                         output_unit, eold, emin, wildcard, gopt_param, SIZE(x0), x0 - xold, xi, conv, &
     271      415038 :                         used_time=t_diff)
     272        2526 :          eold = opt_energy
     273        2526 :          emin = MIN(emin, opt_energy)
     274             : 
     275        2526 :          IF (conv .OR. (its == maxiter)) EXIT
     276             :          ![NB] consistent energies and forces not required for CG, but some line minimizers might set it
     277        2048 :          save_consistent_energy_force = gopt_env%require_consistent_energy_force
     278        2048 :          gopt_env%require_consistent_energy_force = .FALSE.
     279             : 
     280             :          CALL cp_eval_at(gopt_env, x0, opt_energy, xi, gopt_env%force_env%para_env%mepos, &
     281        2048 :                          .FALSE., gopt_env%force_env%para_env)
     282             : 
     283        2048 :          gopt_env%require_consistent_energy_force = save_consistent_energy_force
     284             : 
     285             :          ! Symmetrize coordinates and forces
     286        2048 :          IF (spgr%keep_space_group) THEN
     287           0 :             CALL spgr_apply_rotations_force(spgr, xi)
     288             :          END IF
     289             : 
     290             :          ! Get Conjugate Directions:  updates the searching direction (h)
     291        2048 :          wildcard = "   CG"
     292        2048 :          CALL get_conjugate_direction(gopt_env, Fletcher_Reeves, g, xi, h)
     293             : 
     294             :          ! Symmetrize coordinates and forces
     295        2048 :          IF (spgr%keep_space_group) THEN
     296           0 :             CALL spgr_apply_rotations_force(spgr, g)
     297           0 :             CALL spgr_apply_rotations_force(spgr, h)
     298             :          END IF
     299             : 
     300             :          ! Reset Condition or Steepest Descent Requested
     301             :          ! ABS(DOT_PRODUCT(g, h))/SQRT((DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h))) > res_lim ...
     302             :          IF ((DOT_PRODUCT(g, h)*DOT_PRODUCT(g, h)) > (res_lim*res_lim*DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h)) &
     303      967596 :              .OR. its + 1 <= max_steep_steps) THEN
     304             :             ! Steepest Descent
     305         560 :             wildcard = "   SD"
     306       80728 :             h = -xi
     307             :          END IF
     308      645064 :          g = -xi
     309      648068 :          xi = h
     310             :       END DO
     311             : 
     312         478 :       IF (its == maxiter .AND. (.NOT. conv)) THEN
     313          84 :          CALL print_geo_opt_nc(gopt_env, output_unit)
     314             :       END IF
     315             : 
     316             :       ! Write final particle information and restart, if converged
     317         478 :       IF (PRESENT(do_update)) do_update = conv
     318         478 :       CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
     319             :       CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
     320         478 :                               gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
     321             : 
     322         478 :       DEALLOCATE (xold)
     323         478 :       DEALLOCATE (g)
     324         478 :       DEALLOCATE (h)
     325         478 :       DEALLOCATE (xi)
     326             : 
     327         478 :       CALL timestop(handle)
     328             : 
     329         478 :    END SUBROUTINE cp_cg_main
     330             : 
     331             : END MODULE cg_optimizer

Generated by: LCOV version 1.15