LCOV - code coverage report
Current view: top level - src - qs_ot_eigensolver.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 105 108 97.2 %
Date: 2024-03-28 07:31:50 Functions: 1 1 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 an eigen-space solver for the generalised symmetric eigenvalue problem
      10             : !>      for sparse matrices, needing only multiplications
      11             : !> \author Joost VandeVondele (25.08.2002)
      12             : ! **************************************************************************************************
      13             : MODULE qs_ot_eigensolver
      14             :    USE cp_dbcsr_cholesky,               ONLY: cp_dbcsr_cholesky_decompose,&
      15             :                                               cp_dbcsr_cholesky_invert
      16             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      17             :                                               cp_dbcsr_m_by_n_from_template,&
      18             :                                               cp_fm_to_dbcsr_row_template,&
      19             :                                               dbcsr_copy_columns_hack
      20             :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      21             :                                               cp_fm_type
      22             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      23             :    USE dbcsr_api,                       ONLY: &
      24             :         dbcsr_copy, dbcsr_dot, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, &
      25             :         dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
      26             :    USE kinds,                           ONLY: dp
      27             :    USE preconditioner_types,            ONLY: preconditioner_in_use,&
      28             :                                               preconditioner_type
      29             :    USE qs_mo_methods,                   ONLY: make_basis_sv
      30             :    USE qs_ot,                           ONLY: qs_ot_get_orbitals,&
      31             :                                               qs_ot_get_p,&
      32             :                                               qs_ot_new_preconditioner
      33             :    USE qs_ot_minimizer,                 ONLY: ot_mini
      34             :    USE qs_ot_types,                     ONLY: qs_ot_allocate,&
      35             :                                               qs_ot_destroy,&
      36             :                                               qs_ot_init,&
      37             :                                               qs_ot_settings_init,&
      38             :                                               qs_ot_settings_type,&
      39             :                                               qs_ot_type
      40             : #include "./base/base_uses.f90"
      41             : 
      42             :    IMPLICIT NONE
      43             :    PRIVATE
      44             : 
      45             : ! *** Global parameters ***
      46             : 
      47             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_eigensolver'
      48             : 
      49             : ! *** Public subroutines ***
      50             : 
      51             :    PUBLIC :: ot_eigensolver
      52             : 
      53             : CONTAINS
      54             : 
      55             : ! on input c contains the initial guess (should not be zero !)
      56             : ! on output c spans the subspace
      57             : ! **************************************************************************************************
      58             : !> \brief ...
      59             : !> \param matrix_h ...
      60             : !> \param matrix_s ...
      61             : !> \param matrix_orthogonal_space_fm ...
      62             : !> \param matrix_c_fm ...
      63             : !> \param preconditioner ...
      64             : !> \param eps_gradient ...
      65             : !> \param iter_max ...
      66             : !> \param size_ortho_space ...
      67             : !> \param silent ...
      68             : !> \param ot_settings ...
      69             : ! **************************************************************************************************
      70         836 :    SUBROUTINE ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, &
      71             :                              matrix_c_fm, preconditioner, eps_gradient, &
      72             :                              iter_max, size_ortho_space, silent, ot_settings)
      73             : 
      74             :       TYPE(dbcsr_type), POINTER                          :: matrix_h, matrix_s
      75             :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: matrix_orthogonal_space_fm
      76             :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_c_fm
      77             :       TYPE(preconditioner_type), OPTIONAL, POINTER       :: preconditioner
      78             :       REAL(KIND=dp)                                      :: eps_gradient
      79             :       INTEGER, INTENT(IN)                                :: iter_max
      80             :       INTEGER, INTENT(IN), OPTIONAL                      :: size_ortho_space
      81             :       LOGICAL, INTENT(IN), OPTIONAL                      :: silent
      82             :       TYPE(qs_ot_settings_type), INTENT(IN), OPTIONAL    :: ot_settings
      83             : 
      84             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ot_eigensolver'
      85             :       INTEGER, PARAMETER                                 :: max_iter_inner_loop = 40
      86             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
      87             : 
      88             :       INTEGER                                            :: handle, ieigensolver, iter_total, k, n, &
      89             :                                                             ortho_k, ortho_space_k, output_unit
      90             :       LOGICAL                                            :: energy_only, my_silent, ortho
      91             :       REAL(KIND=dp)                                      :: delta, energy
      92         418 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_hc
      93             :       TYPE(dbcsr_type), POINTER                          :: matrix_buf1_ortho, matrix_buf2_ortho, &
      94             :                                                             matrix_c, matrix_orthogonal_space, &
      95             :                                                             matrix_os_ortho, matrix_s_ortho
      96         418 :       TYPE(qs_ot_type), DIMENSION(:), POINTER            :: qs_ot_env
      97             : 
      98         418 :       CALL timeset(routineN, handle)
      99             : 
     100         418 :       output_unit = cp_logger_get_default_io_unit()
     101             : 
     102         418 :       IF (PRESENT(silent)) THEN
     103         116 :          my_silent = silent
     104             :       ELSE
     105             :          my_silent = .FALSE.
     106             :       END IF
     107             : 
     108         418 :       NULLIFY (matrix_c) ! fm->dbcsr
     109             : 
     110         418 :       CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k) ! fm->dbcsr
     111         418 :       ALLOCATE (matrix_c)
     112         418 :       CALL cp_fm_to_dbcsr_row_template(matrix_c, fm_in=matrix_c_fm, template=matrix_h)
     113             : 
     114         418 :       iter_total = 0
     115             : 
     116             :       outer_scf: DO
     117             : 
     118             :          NULLIFY (qs_ot_env)
     119             : 
     120         654 :          NULLIFY (matrix_s_ortho)
     121         654 :          NULLIFY (matrix_os_ortho)
     122         654 :          NULLIFY (matrix_buf1_ortho)
     123         654 :          NULLIFY (matrix_buf2_ortho)
     124         654 :          NULLIFY (matrix_orthogonal_space)
     125             : 
     126         654 :          ALLOCATE (qs_ot_env(1))
     127        1308 :          ALLOCATE (matrix_hc(1))
     128         654 :          NULLIFY (matrix_hc(1)%matrix)
     129         654 :          CALL dbcsr_init_p(matrix_hc(1)%matrix)
     130         654 :          CALL dbcsr_copy(matrix_hc(1)%matrix, matrix_c, 'matrix_hc')
     131             : 
     132         654 :          ortho = .FALSE.
     133         654 :          IF (PRESENT(matrix_orthogonal_space_fm)) ortho = .TRUE.
     134             : 
     135             :          ! decide settings
     136         654 :          IF (PRESENT(ot_settings)) THEN
     137         118 :             qs_ot_env(1)%settings = ot_settings
     138             :          ELSE
     139         536 :             CALL qs_ot_settings_init(qs_ot_env(1)%settings)
     140             :             ! overwrite defaults
     141         536 :             qs_ot_env(1)%settings%ds_min = 0.10_dp
     142             :          END IF
     143             : 
     144         654 :          IF (ortho) THEN
     145         536 :             ALLOCATE (matrix_orthogonal_space)
     146         536 :             CALL cp_fm_to_dbcsr_row_template(matrix_orthogonal_space, fm_in=matrix_orthogonal_space_fm, template=matrix_h)
     147         536 :             CALL cp_fm_get_info(matrix_orthogonal_space_fm, ncol_global=ortho_space_k)
     148             : 
     149         536 :             IF (PRESENT(size_ortho_space)) ortho_space_k = size_ortho_space
     150         536 :             ortho_k = ortho_space_k + k
     151             :          ELSE
     152         118 :             ortho_k = k
     153             :          END IF
     154             : 
     155             :          ! allocate
     156         654 :          CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct, ortho_k=ortho_k)
     157             : 
     158         654 :          IF (ortho) THEN
     159             :             ! construct an initial guess that is orthogonal to matrix_orthogonal_space
     160             : 
     161         536 :             CALL dbcsr_init_p(matrix_s_ortho)
     162         536 :             CALL dbcsr_copy(matrix_s_ortho, matrix_orthogonal_space, name="matrix_s_ortho")
     163             : 
     164         536 :             CALL dbcsr_init_p(matrix_os_ortho)
     165             :             CALL cp_dbcsr_m_by_n_from_template(matrix_os_ortho, template=matrix_h, m=ortho_space_k, n=ortho_space_k, &
     166         536 :                                                sym=dbcsr_type_no_symmetry)
     167             : 
     168         536 :             CALL dbcsr_init_p(matrix_buf1_ortho)
     169             :             CALL cp_dbcsr_m_by_n_from_template(matrix_buf1_ortho, template=matrix_h, m=ortho_space_k, n=k, &
     170         536 :                                                sym=dbcsr_type_no_symmetry)
     171             : 
     172         536 :             CALL dbcsr_init_p(matrix_buf2_ortho)
     173             :             CALL cp_dbcsr_m_by_n_from_template(matrix_buf2_ortho, template=matrix_h, m=ortho_space_k, n=k, &
     174         536 :                                                sym=dbcsr_type_no_symmetry)
     175             : 
     176             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_orthogonal_space, &
     177         536 :                                 0.0_dp, matrix_s_ortho)
     178             :             CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_s_ortho, &
     179         536 :                                 rzero, matrix_os_ortho)
     180             : 
     181             :             CALL cp_dbcsr_cholesky_decompose(matrix_os_ortho, &
     182         536 :                                              para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
     183             :             CALL cp_dbcsr_cholesky_invert(matrix_os_ortho, &
     184             :                                           para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env, &
     185         536 :                                           upper_to_full=.TRUE.)
     186             : 
     187             :             CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_c, &
     188         536 :                                 rzero, matrix_buf1_ortho)
     189             :             CALL dbcsr_multiply('N', 'N', rone, matrix_os_ortho, matrix_buf1_ortho, &
     190         536 :                                 rzero, matrix_buf2_ortho)
     191             :             CALL dbcsr_multiply('N', 'N', -rone, matrix_s_ortho, matrix_buf2_ortho, &
     192         536 :                                 rone, matrix_c)
     193             : 
     194             :             ! make matrix_c0 an orthogonal basis, matrix_c contains sc0
     195         536 :             CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
     196             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
     197         536 :                                 0.0_dp, matrix_c)
     198             : 
     199             :             CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, matrix_c, &
     200         536 :                                qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
     201             : 
     202             :             ! copy sc0 and matrix_s_ortho in qs_ot_env(1)%matrix_sc0
     203             :             !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_s_ortho,ortho_space_k,1,1)
     204             :             CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_s_ortho, ortho_space_k, 1, 1, &
     205         536 :                                          para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
     206             :             !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_c,k,1,ortho_space_k+1)
     207             :             CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_c, k, 1, ortho_space_k + 1, &
     208         536 :                                          para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
     209             : 
     210         536 :             CALL dbcsr_release_p(matrix_buf1_ortho)
     211         536 :             CALL dbcsr_release_p(matrix_buf2_ortho)
     212         536 :             CALL dbcsr_release_p(matrix_os_ortho)
     213         536 :             CALL dbcsr_release_p(matrix_s_ortho)
     214             : 
     215             :          ELSE
     216             : 
     217             :             ! set c0,sc0
     218         118 :             CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
     219             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
     220         118 :                                 0.0_dp, qs_ot_env(1)%matrix_sc0)
     221             : 
     222             :             CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, qs_ot_env(1)%matrix_sc0, &
     223         118 :                                qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
     224             :          END IF
     225             : 
     226             :          ! init
     227         654 :          CALL qs_ot_init(qs_ot_env(1))
     228         654 :          energy_only = qs_ot_env(1)%energy_only
     229             : 
     230             :          ! set x
     231         654 :          CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
     232         654 :          CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
     233             : 
     234             :          ! get c
     235         654 :          CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
     236         654 :          CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
     237             : 
     238             :          ! if present preconditioner, use it
     239             : 
     240         654 :          IF (PRESENT(preconditioner)) THEN
     241         654 :             IF (ASSOCIATED(preconditioner)) THEN
     242         294 :                IF (preconditioner_in_use(preconditioner)) THEN
     243         294 :                   CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
     244             :                ELSE
     245             :                   ! we should presumably make one
     246             :                END IF
     247             :             END IF
     248             :          END IF
     249             : 
     250             :          ! *** Eigensolver loop ***
     251             :          ieigensolver = 0
     252             :          eigensolver_loop: DO
     253             : 
     254       14676 :             ieigensolver = ieigensolver + 1
     255       14676 :             iter_total = iter_total + 1
     256             : 
     257             :             ! the energy is cHc, the gradient is 2*H*c
     258             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c, &
     259       14676 :                                 0.0_dp, matrix_hc(1)%matrix)
     260       14676 :             CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy)
     261       14676 :             IF (.NOT. energy_only) THEN
     262        7638 :                CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
     263             :             END IF
     264             : 
     265       14676 :             qs_ot_env(1)%etotal = energy
     266       14676 :             CALL ot_mini(qs_ot_env, matrix_hc)
     267       14676 :             delta = qs_ot_env(1)%delta
     268       14676 :             energy_only = qs_ot_env(1)%energy_only
     269             : 
     270             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_x, &
     271       14676 :                                 0.0_dp, qs_ot_env(1)%matrix_sx)
     272             : 
     273       14676 :             CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
     274       14676 :             CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
     275             : 
     276             :             ! exit on convergence or if maximum of inner loop  cycles is reached
     277       14676 :             IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
     278             :             ! exit if total number of steps is reached, but not during a line search step
     279       14676 :             IF (iter_total >= iter_max .AND. qs_ot_env(1)%OT_METHOD_FULL /= "OT LS") EXIT eigensolver_loop
     280             : 
     281             :          END DO eigensolver_loop
     282             : 
     283         654 :          CALL qs_ot_destroy(qs_ot_env(1))
     284         654 :          DEALLOCATE (qs_ot_env)
     285         654 :          CALL dbcsr_release_p(matrix_hc(1)%matrix)
     286         654 :          DEALLOCATE (matrix_hc)
     287         654 :          CALL dbcsr_release_p(matrix_orthogonal_space)
     288             : 
     289         654 :          IF (delta < eps_gradient) THEN
     290         368 :             IF ((output_unit > 0) .AND. .NOT. my_silent) THEN
     291             :                WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
     292         155 :                   "OT| Eigensolver reached convergence in ", iter_total, " iterations"
     293             :             END IF
     294             :             EXIT outer_scf
     295             :          END IF
     296         286 :          IF (iter_total >= iter_max) THEN
     297          50 :             IF (output_unit > 0) THEN
     298          25 :                IF (my_silent) THEN
     299          25 :                   WRITE (output_unit, "(A,T60,E20.10)") "  WARNING OT eigensolver did not converge: current gradient", delta
     300             :                ELSE
     301           0 :                   WRITE (output_unit, *) "WARNING : did not converge in ot_eigensolver"
     302           0 :                   WRITE (output_unit, *) "number of iterations ", iter_total, " exceeded maximum"
     303           0 :                   WRITE (output_unit, *) "current gradient / target gradient", delta, " / ", eps_gradient
     304             :                END IF
     305             :             END IF
     306             :             EXIT outer_scf
     307             :          END IF
     308             : 
     309             :       END DO outer_scf
     310             : 
     311         418 :       CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm) ! fm->dbcsr
     312         418 :       CALL dbcsr_release_p(matrix_c) ! fm->dbcsr
     313             : 
     314         418 :       CALL timestop(handle)
     315             : 
     316         418 :    END SUBROUTINE ot_eigensolver
     317             : 
     318             : END MODULE qs_ot_eigensolver

Generated by: LCOV version 1.15