LCOV - code coverage report
Current view: top level - src - qs_ot_eigensolver.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.2 % 109 106
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 1 1

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

Generated by: LCOV version 2.0-1