LCOV - code coverage report
Current view: top level - src - preconditioner_makes.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 91.1 % 305 278
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 8 8

            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 computes preconditioners, and implements methods to apply them
      10              : !>      currently used in qs_ot
      11              : !> \par History
      12              : !>      - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse)
      13              : !> \author Joost VandeVondele (09.2002)
      14              : ! **************************************************************************************************
      15              : MODULE preconditioner_makes
      16              :    USE arnoldi_api,                     ONLY: arnoldi_env_type,&
      17              :                                               arnoldi_ev,&
      18              :                                               deallocate_arnoldi_env,&
      19              :                                               get_selected_ritz_val,&
      20              :                                               get_selected_ritz_vector,&
      21              :                                               set_arnoldi_initial_vector,&
      22              :                                               setup_arnoldi_env
      23              :    USE cp_dbcsr_api,                    ONLY: &
      24              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_p_type, &
      25              :         dbcsr_release, dbcsr_type, dbcsr_type_symmetric
      26              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_add_on_diag
      27              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      28              :                                               cp_dbcsr_m_by_n_from_template,&
      29              :                                               cp_dbcsr_sm_fm_multiply,&
      30              :                                               cp_fm_to_dbcsr_row_template
      31              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      32              :                                               cp_fm_triangular_invert,&
      33              :                                               cp_fm_triangular_multiply,&
      34              :                                               cp_fm_uplo_to_full
      35              :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose,&
      36              :                                               cp_fm_cholesky_reduce,&
      37              :                                               cp_fm_cholesky_restore
      38              :    USE cp_fm_diag,                      ONLY: choose_eigv_solver
      39              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      40              :                                               cp_fm_struct_release,&
      41              :                                               cp_fm_struct_type
      42              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      43              :                                               cp_fm_get_diag,&
      44              :                                               cp_fm_get_info,&
      45              :                                               cp_fm_release,&
      46              :                                               cp_fm_to_fm,&
      47              :                                               cp_fm_type
      48              :    USE input_constants,                 ONLY: &
      49              :         cholesky_inverse, cholesky_reduce, ot_precond_full_all, ot_precond_full_kinetic, &
      50              :         ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_s_inverse, &
      51              :         ot_precond_solver_default, ot_precond_solver_inv_chol
      52              :    USE kinds,                           ONLY: dp
      53              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      54              :    USE preconditioner_types,            ONLY: preconditioner_type
      55              : #include "./base/base_uses.f90"
      56              : 
      57              :    IMPLICIT NONE
      58              : 
      59              :    PRIVATE
      60              : 
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_makes'
      62              : 
      63              :    PUBLIC :: make_preconditioner_matrix
      64              : 
      65              : CONTAINS
      66              : 
      67              : ! **************************************************************************************************
      68              : !> \brief ...
      69              : !> \param preconditioner_env ...
      70              : !> \param matrix_h ...
      71              : !> \param matrix_s ...
      72              : !> \param matrix_t ...
      73              : !> \param mo_coeff ...
      74              : !> \param energy_homo ...
      75              : !> \param eigenvalues_ot ...
      76              : !> \param energy_gap ...
      77              : !> \param my_solver_type ...
      78              : ! **************************************************************************************************
      79         8500 :    SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
      80         8500 :                                          energy_homo, eigenvalues_ot, energy_gap, &
      81              :                                          my_solver_type)
      82              :       TYPE(preconditioner_type)                          :: preconditioner_env
      83              :       TYPE(dbcsr_type), POINTER                          :: matrix_h
      84              :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s, matrix_t
      85              :       TYPE(cp_fm_type), INTENT(IN)                       :: mo_coeff
      86              :       REAL(KIND=dp)                                      :: energy_homo
      87              :       REAL(KIND=dp), DIMENSION(:)                        :: eigenvalues_ot
      88              :       REAL(KIND=dp)                                      :: energy_gap
      89              :       INTEGER                                            :: my_solver_type
      90              : 
      91              :       INTEGER                                            :: precon_type
      92              : 
      93         8500 :       precon_type = preconditioner_env%in_use
      94           38 :       SELECT CASE (precon_type)
      95              :       CASE (ot_precond_full_single)
      96           38 :          IF (my_solver_type .NE. ot_precond_solver_default) &
      97            0 :             CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
      98           38 :          IF (PRESENT(matrix_s)) THEN
      99              :             CALL make_full_single(preconditioner_env, preconditioner_env%fm, &
     100           32 :                                   matrix_h, matrix_s, energy_homo, energy_gap)
     101              :          ELSE
     102              :             CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm, &
     103            6 :                                         matrix_h, energy_homo, energy_gap)
     104              :          END IF
     105              : 
     106              :       CASE (ot_precond_s_inverse)
     107           72 :          IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
     108           72 :          IF (.NOT. PRESENT(matrix_s)) &
     109            0 :             CPABORT("Type for S=1 not implemented")
     110           72 :          CALL make_full_s_inverse(preconditioner_env, matrix_s)
     111              : 
     112              :       CASE (ot_precond_full_kinetic)
     113         1305 :          IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
     114         1305 :          IF (.NOT. (PRESENT(matrix_s) .AND. PRESENT(matrix_t))) &
     115            0 :             CPABORT("Type for S=1 not implemented")
     116         1305 :          CALL make_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
     117              :       CASE (ot_precond_full_single_inverse)
     118         4045 :          IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
     119              :          CALL make_full_single_inverse(preconditioner_env, mo_coeff, matrix_h, energy_gap, &
     120         4045 :                                        matrix_s=matrix_s)
     121              :       CASE (ot_precond_full_all)
     122         3040 :          IF (my_solver_type .NE. ot_precond_solver_default) THEN
     123            0 :             CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
     124              :          END IF
     125         3040 :          IF (PRESENT(matrix_s)) THEN
     126              :             CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
     127         2958 :                                eigenvalues_ot, energy_gap)
     128              :          ELSE
     129              :             CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
     130           82 :                                      eigenvalues_ot, energy_gap)
     131              :          END IF
     132              : 
     133              :       CASE DEFAULT
     134         8500 :          CPABORT("Type not implemented")
     135              :       END SELECT
     136              : 
     137         8500 :    END SUBROUTINE make_preconditioner_matrix
     138              : 
     139              : ! **************************************************************************************************
     140              : !> \brief Simply takes the overlap matrix as preconditioner
     141              : !> \param preconditioner_env ...
     142              : !> \param matrix_s ...
     143              : ! **************************************************************************************************
     144           72 :    SUBROUTINE make_full_s_inverse(preconditioner_env, matrix_s)
     145              :       TYPE(preconditioner_type)                          :: preconditioner_env
     146              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     147              : 
     148              :       CHARACTER(len=*), PARAMETER :: routineN = 'make_full_s_inverse'
     149              : 
     150              :       INTEGER                                            :: handle
     151              : 
     152           72 :       CALL timeset(routineN, handle)
     153              : 
     154           72 :       CPASSERT(ASSOCIATED(matrix_s))
     155              : 
     156           72 :       IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
     157           72 :          ALLOCATE (preconditioner_env%sparse_matrix)
     158              :       END IF
     159           72 :       CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic")
     160              : 
     161           72 :       CALL timestop(handle)
     162              : 
     163           72 :    END SUBROUTINE make_full_s_inverse
     164              : 
     165              : ! **************************************************************************************************
     166              : !> \brief kinetic matrix+shift*overlap as preconditioner. Cheap but could
     167              : !>        be better
     168              : !> \param preconditioner_env ...
     169              : !> \param matrix_t ...
     170              : !> \param matrix_s ...
     171              : !> \param energy_gap ...
     172              : ! **************************************************************************************************
     173         1305 :    SUBROUTINE make_full_kinetic(preconditioner_env, matrix_t, matrix_s, &
     174              :                                 energy_gap)
     175              :       TYPE(preconditioner_type)                          :: preconditioner_env
     176              :       TYPE(dbcsr_type), POINTER                          :: matrix_t, matrix_s
     177              :       REAL(KIND=dp)                                      :: energy_gap
     178              : 
     179              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'make_full_kinetic'
     180              : 
     181              :       INTEGER                                            :: handle
     182              :       REAL(KIND=dp)                                      :: shift
     183              : 
     184         1305 :       CALL timeset(routineN, handle)
     185              : 
     186         1305 :       CPASSERT(ASSOCIATED(matrix_t))
     187         1305 :       CPASSERT(ASSOCIATED(matrix_s))
     188              : 
     189         1305 :       IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
     190         1303 :          ALLOCATE (preconditioner_env%sparse_matrix)
     191              :       END IF
     192         1305 :       CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic")
     193              : 
     194         1305 :       shift = MAX(0.0_dp, energy_gap)
     195              : 
     196              :       CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, &
     197         1305 :                      alpha_scalar=1.0_dp, beta_scalar=shift)
     198              : 
     199         1305 :       CALL timestop(handle)
     200              : 
     201         1305 :    END SUBROUTINE make_full_kinetic
     202              : 
     203              : ! **************************************************************************************************
     204              : !> \brief full_single_preconditioner
     205              : !> \param preconditioner_env ...
     206              : !> \param fm ...
     207              : !> \param matrix_h ...
     208              : !> \param matrix_s ...
     209              : !> \param energy_homo ...
     210              : !> \param energy_gap ...
     211              : ! **************************************************************************************************
     212           32 :    SUBROUTINE make_full_single(preconditioner_env, fm, matrix_h, matrix_s, &
     213              :                                energy_homo, energy_gap)
     214              :       TYPE(preconditioner_type)                          :: preconditioner_env
     215              :       TYPE(cp_fm_type), POINTER                          :: fm
     216              :       TYPE(dbcsr_type), POINTER                          :: matrix_h, matrix_s
     217              :       REAL(KIND=dp)                                      :: energy_homo, energy_gap
     218              : 
     219              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'make_full_single'
     220              : 
     221              :       INTEGER                                            :: handle, i, n
     222           32 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: evals
     223              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     224              :       TYPE(cp_fm_type)                                   :: fm_h, fm_s
     225              : 
     226           32 :       CALL timeset(routineN, handle)
     227              : 
     228           32 :       NULLIFY (fm_struct_tmp, evals)
     229              : 
     230           32 :       IF (ASSOCIATED(fm)) THEN
     231            0 :          CALL cp_fm_release(fm)
     232            0 :          DEALLOCATE (fm)
     233              :          NULLIFY (fm)
     234              :       END IF
     235           32 :       CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
     236           96 :       ALLOCATE (evals(n))
     237              : 
     238              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
     239              :                                context=preconditioner_env%ctxt, &
     240           32 :                                para_env=preconditioner_env%para_env)
     241           32 :       ALLOCATE (fm)
     242           32 :       CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
     243           32 :       CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
     244           32 :       CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
     245           32 :       CALL cp_fm_struct_release(fm_struct_tmp)
     246              : 
     247           32 :       CALL copy_dbcsr_to_fm(matrix_h, fm_h)
     248           32 :       CALL copy_dbcsr_to_fm(matrix_s, fm_s)
     249           32 :       CALL cp_fm_cholesky_decompose(fm_s)
     250              : 
     251           32 :       SELECT CASE (preconditioner_env%cholesky_use)
     252              :       CASE (cholesky_inverse)
     253              : ! if cho inverse
     254            0 :          CALL cp_fm_triangular_invert(fm_s)
     255            0 :          CALL cp_fm_uplo_to_full(fm_h, fm)
     256              : 
     257              :          CALL cp_fm_triangular_multiply(fm_s, fm_h, side="R", transpose_tr=.FALSE., &
     258            0 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
     259              :          CALL cp_fm_triangular_multiply(fm_s, fm_h, side="L", transpose_tr=.TRUE., &
     260            0 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
     261              :       CASE (cholesky_reduce)
     262           32 :          CALL cp_fm_cholesky_reduce(fm_h, fm_s)
     263              :       CASE DEFAULT
     264           32 :          CPABORT("cholesky type not implemented")
     265              :       END SELECT
     266              : 
     267           32 :       CALL choose_eigv_solver(fm_h, fm, evals)
     268              : 
     269           32 :       SELECT CASE (preconditioner_env%cholesky_use)
     270              :       CASE (cholesky_inverse)
     271              :          CALL cp_fm_triangular_multiply(fm_s, fm, side="L", transpose_tr=.FALSE., &
     272            0 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
     273            0 :          DO i = 1, n
     274            0 :             evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
     275              :          END DO
     276            0 :          CALL cp_fm_to_fm(fm, fm_h)
     277              :       CASE (cholesky_reduce)
     278           32 :          CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE")
     279          568 :          DO i = 1, n
     280          568 :             evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
     281              :          END DO
     282           64 :          CALL cp_fm_to_fm(fm_h, fm)
     283              :       END SELECT
     284              : 
     285           32 :       CALL cp_fm_column_scale(fm, evals)
     286           32 :       CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
     287           32 :       CALL cp_fm_to_fm(fm_s, fm)
     288              : 
     289           32 :       DEALLOCATE (evals)
     290           32 :       CALL cp_fm_release(fm_h)
     291           32 :       CALL cp_fm_release(fm_s)
     292              : 
     293           32 :       CALL timestop(handle)
     294              : 
     295           96 :    END SUBROUTINE make_full_single
     296              : 
     297              : ! **************************************************************************************************
     298              : !> \brief full single in the orthonormal basis
     299              : !> \param preconditioner_env ...
     300              : !> \param fm ...
     301              : !> \param matrix_h ...
     302              : !> \param energy_homo ...
     303              : !> \param energy_gap ...
     304              : ! **************************************************************************************************
     305            6 :    SUBROUTINE make_full_single_ortho(preconditioner_env, fm, matrix_h, &
     306              :                                      energy_homo, energy_gap)
     307              :       TYPE(preconditioner_type)                          :: preconditioner_env
     308              :       TYPE(cp_fm_type), POINTER                          :: fm
     309              :       TYPE(dbcsr_type), POINTER                          :: matrix_h
     310              :       REAL(KIND=dp)                                      :: energy_homo, energy_gap
     311              : 
     312              :       CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_ortho'
     313              : 
     314              :       INTEGER                                            :: handle, i, n
     315            6 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: evals
     316              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     317              :       TYPE(cp_fm_type)                                   :: fm_h, fm_s
     318              : 
     319            6 :       CALL timeset(routineN, handle)
     320            6 :       NULLIFY (fm_struct_tmp, evals)
     321              : 
     322            6 :       IF (ASSOCIATED(fm)) THEN
     323            0 :          CALL cp_fm_release(fm)
     324            0 :          DEALLOCATE (fm)
     325              :          NULLIFY (fm)
     326              :       END IF
     327            6 :       CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
     328           18 :       ALLOCATE (evals(n))
     329              : 
     330              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
     331              :                                context=preconditioner_env%ctxt, &
     332            6 :                                para_env=preconditioner_env%para_env)
     333            6 :       ALLOCATE (fm)
     334            6 :       CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
     335            6 :       CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
     336            6 :       CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
     337            6 :       CALL cp_fm_struct_release(fm_struct_tmp)
     338              : 
     339            6 :       CALL copy_dbcsr_to_fm(matrix_h, fm_h)
     340              : 
     341            6 :       CALL choose_eigv_solver(fm_h, fm, evals)
     342          282 :       DO i = 1, n
     343          282 :          evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
     344              :       END DO
     345            6 :       CALL cp_fm_to_fm(fm, fm_h)
     346            6 :       CALL cp_fm_column_scale(fm, evals)
     347            6 :       CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
     348            6 :       CALL cp_fm_to_fm(fm_s, fm)
     349              : 
     350            6 :       DEALLOCATE (evals)
     351            6 :       CALL cp_fm_release(fm_h)
     352            6 :       CALL cp_fm_release(fm_s)
     353              : 
     354            6 :       CALL timestop(handle)
     355              : 
     356           18 :    END SUBROUTINE make_full_single_ortho
     357              : 
     358              : ! **************************************************************************************************
     359              : !> \brief generates a state by state preconditioner based on the full hamiltonian matrix
     360              : !> \param preconditioner_env ...
     361              : !> \param matrix_c0 ...
     362              : !> \param matrix_h ...
     363              : !> \param matrix_s ...
     364              : !> \param c0_evals ...
     365              : !> \param energy_gap should be a slight underestimate of the physical energy gap for almost all systems
     366              : !>      the c0 are already ritz states of (h,s)
     367              : !> \par History
     368              : !>      10.2006 made more stable [Joost VandeVondele]
     369              : !> \note
     370              : !>      includes error estimate on the hamiltonian matrix to result in a stable preconditioner
     371              : !>      a preconditioner for each eigenstate i is generated by keeping the factorized form
     372              : !>      U diag( something i ) U^T. It is important to only precondition in the subspace orthogonal to c0.
     373              : !>      not only is it the only part that matters, it also simplifies the computation of
     374              : !>      the lagrangian multipliers in the OT minimization  (i.e. if the c0 here is different
     375              : !>      from the c0 used in the OT setup, there will be a bug).
     376              : ! **************************************************************************************************
     377         2958 :    SUBROUTINE make_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap)
     378              :       TYPE(preconditioner_type)                          :: preconditioner_env
     379              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_c0
     380              :       TYPE(dbcsr_type), POINTER                          :: matrix_h, matrix_s
     381              :       REAL(KIND=dp), DIMENSION(:)                        :: c0_evals
     382              :       REAL(KIND=dp)                                      :: energy_gap
     383              : 
     384              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'make_full_all'
     385              :       REAL(KIND=dp), PARAMETER                           :: fudge_factor = 0.25_dp, &
     386              :                                                             lambda_base = 10.0_dp
     387              : 
     388              :       INTEGER                                            :: handle, k, n
     389              :       REAL(KIND=dp)                                      :: error_estimate, lambda
     390         2958 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: diag, norms, shifted_evals
     391              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     392              :       TYPE(cp_fm_type)                                   :: matrix_hc0, matrix_left, matrix_s1, &
     393              :                                                             matrix_s2, matrix_sc0, matrix_shc0, &
     394              :                                                             matrix_tmp, ortho
     395              :       TYPE(cp_fm_type), POINTER                          :: matrix_pre
     396              : 
     397         2958 :       CALL timeset(routineN, handle)
     398              : 
     399         2958 :       IF (ASSOCIATED(preconditioner_env%fm)) THEN
     400            0 :          CALL cp_fm_release(preconditioner_env%fm)
     401            0 :          DEALLOCATE (preconditioner_env%fm)
     402              :          NULLIFY (preconditioner_env%fm)
     403              :       END IF
     404         2958 :       CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
     405              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
     406              :                                context=preconditioner_env%ctxt, &
     407         2958 :                                para_env=preconditioner_env%para_env)
     408         2958 :       ALLOCATE (preconditioner_env%fm)
     409         2958 :       CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
     410         2958 :       CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho")
     411         2958 :       CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
     412         2958 :       CALL cp_fm_struct_release(fm_struct_tmp)
     413         8874 :       ALLOCATE (preconditioner_env%full_evals(n))
     414         8768 :       ALLOCATE (preconditioner_env%occ_evals(k))
     415              : 
     416              :       ! 0) cholesky decompose the overlap matrix, if this fails the basis is singular,
     417              :       !    more than EPS_DEFAULT
     418         2958 :       CALL copy_dbcsr_to_fm(matrix_s, ortho)
     419         2958 :       CALL cp_fm_cholesky_decompose(ortho)
     420              : ! if cho inverse
     421         2958 :       IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN
     422            0 :          CALL cp_fm_triangular_invert(ortho)
     423              :       END IF
     424              :       ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
     425              :       !    possibly shifted by an amount lambda,
     426              :       !    and the same spectrum as the original H matrix in the space orthogonal to the C0
     427              :       !    with P=C0 C0 ^ T
     428              :       !    (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
     429              :       !    we exploit that the C0 are already the ritz states of H
     430         2958 :       CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
     431         2958 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k)
     432         2958 :       CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
     433         2958 :       CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
     434              : 
     435              :       ! An aside, try to estimate the error on the ritz values, we'll need it later on
     436         2958 :       CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0")
     437              : 
     438         2958 :       SELECT CASE (preconditioner_env%cholesky_use)
     439              :       CASE (cholesky_inverse)
     440              : ! if cho inverse
     441            0 :          CALL cp_fm_to_fm(matrix_hc0, matrix_shc0)
     442              :          CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.TRUE., &
     443            0 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp)
     444              :       CASE (cholesky_reduce)
     445         2958 :          CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T")
     446              :       CASE DEFAULT
     447         2958 :          CPABORT("cholesky type not implemented")
     448              :       END SELECT
     449              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
     450              :                                context=preconditioner_env%ctxt, &
     451         2958 :                                para_env=preconditioner_env%para_env)
     452         2958 :       CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
     453         2958 :       CALL cp_fm_struct_release(fm_struct_tmp)
     454              :       ! since we only use diagonal elements this is a bit of a waste
     455         2958 :       CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1)
     456         5810 :       ALLOCATE (diag(k))
     457         2958 :       CALL cp_fm_get_diag(matrix_s1, diag)
     458        19818 :       error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
     459         2958 :       DEALLOCATE (diag)
     460         2958 :       CALL cp_fm_release(matrix_s1)
     461         2958 :       CALL cp_fm_release(matrix_shc0)
     462              :       ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
     463              :       ! is small enough. A large error combined with a small energy gap would otherwise lead to
     464              :       ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
     465              :       ! aggressively
     466         2958 :       preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
     467         2958 :       CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
     468         2958 :       matrix_pre => preconditioner_env%fm
     469         2958 :       CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
     470              :       ! tmp = H ( 1 - PS )
     471         2958 :       CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
     472              : 
     473              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
     474              :                                context=preconditioner_env%ctxt, &
     475         2958 :                                para_env=preconditioner_env%para_env)
     476         2958 :       CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
     477         2958 :       CALL cp_fm_struct_release(fm_struct_tmp)
     478         2958 :       CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
     479              :       ! tmp = (1 - PS)^T H (1-PS)
     480         2958 :       CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
     481         2958 :       CALL cp_fm_release(matrix_left)
     482              : 
     483         5810 :       ALLOCATE (shifted_evals(k))
     484         2958 :       lambda = lambda_base + error_estimate
     485        16860 :       shifted_evals = c0_evals - lambda
     486         2958 :       CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
     487         2958 :       CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
     488         2958 :       CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
     489              : 
     490              :       ! 2) diagonalize this operator
     491         2958 :       SELECT CASE (preconditioner_env%cholesky_use)
     492              :       CASE (cholesky_inverse)
     493              :          CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="R", transpose_tr=.FALSE., &
     494            0 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
     495              :          CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="L", transpose_tr=.TRUE., &
     496            0 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
     497              :       CASE (cholesky_reduce)
     498         2958 :          CALL cp_fm_cholesky_reduce(matrix_tmp, ortho)
     499              :       END SELECT
     500         2958 :       CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
     501         2958 :       SELECT CASE (preconditioner_env%cholesky_use)
     502              :       CASE (cholesky_inverse)
     503              :          CALL cp_fm_triangular_multiply(ortho, matrix_pre, side="L", transpose_tr=.FALSE., &
     504            0 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
     505            0 :          CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
     506              :       CASE (cholesky_reduce)
     507         2958 :          CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE")
     508         5916 :          CALL cp_fm_to_fm(matrix_tmp, matrix_pre)
     509              :       END SELECT
     510              : 
     511              :       ! test that the subspace remained conserved
     512              :       IF (.FALSE.) THEN
     513              :          CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
     514              :                                   context=preconditioner_env%ctxt, &
     515              :                                   para_env=preconditioner_env%para_env)
     516              :          CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
     517              :          CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
     518              :          CALL cp_fm_struct_release(fm_struct_tmp)
     519              :          ALLOCATE (norms(k))
     520              :          CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
     521              :          CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
     522              :          WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
     523              :          DEALLOCATE (norms)
     524              :          CALL cp_fm_release(matrix_s1)
     525              :          CALL cp_fm_release(matrix_s2)
     526              :       END IF
     527              : 
     528              :       ! 3) replace the lowest k evals and evecs with what they should be
     529        16860 :       preconditioner_env%occ_evals = c0_evals
     530              :       ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all)
     531        16860 :       preconditioner_env%full_evals(1:k) = c0_evals
     532         2958 :       CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
     533              : 
     534         2958 :       CALL cp_fm_release(matrix_sc0)
     535         2958 :       CALL cp_fm_release(matrix_hc0)
     536         2958 :       CALL cp_fm_release(ortho)
     537         2958 :       CALL cp_fm_release(matrix_tmp)
     538         2958 :       DEALLOCATE (shifted_evals)
     539         2958 :       CALL timestop(handle)
     540              : 
     541        23664 :    END SUBROUTINE make_full_all
     542              : 
     543              : ! **************************************************************************************************
     544              : !> \brief full all in the orthonormal basis
     545              : !> \param preconditioner_env ...
     546              : !> \param matrix_c0 ...
     547              : !> \param matrix_h ...
     548              : !> \param c0_evals ...
     549              : !> \param energy_gap ...
     550              : ! **************************************************************************************************
     551           82 :    SUBROUTINE make_full_all_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals, energy_gap)
     552              : 
     553              :       TYPE(preconditioner_type)                          :: preconditioner_env
     554              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_c0
     555              :       TYPE(dbcsr_type), POINTER                          :: matrix_h
     556              :       REAL(KIND=dp), DIMENSION(:)                        :: c0_evals
     557              :       REAL(KIND=dp)                                      :: energy_gap
     558              : 
     559              :       CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho'
     560              :       REAL(KIND=dp), PARAMETER                           :: fudge_factor = 0.25_dp, &
     561              :                                                             lambda_base = 10.0_dp
     562              : 
     563              :       INTEGER                                            :: handle, k, n
     564              :       REAL(KIND=dp)                                      :: error_estimate, lambda
     565           82 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: diag, norms, shifted_evals
     566              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     567              :       TYPE(cp_fm_type)                                   :: matrix_hc0, matrix_left, matrix_s1, &
     568              :                                                             matrix_s2, matrix_sc0, matrix_tmp
     569              :       TYPE(cp_fm_type), POINTER                          :: matrix_pre
     570              : 
     571           82 :       CALL timeset(routineN, handle)
     572              : 
     573           82 :       IF (ASSOCIATED(preconditioner_env%fm)) THEN
     574            0 :          CALL cp_fm_release(preconditioner_env%fm)
     575            0 :          DEALLOCATE (preconditioner_env%fm)
     576              :          NULLIFY (preconditioner_env%fm)
     577              :       END IF
     578           82 :       CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
     579              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
     580              :                                context=preconditioner_env%ctxt, &
     581           82 :                                para_env=preconditioner_env%para_env)
     582           82 :       ALLOCATE (preconditioner_env%fm)
     583           82 :       CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
     584           82 :       CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
     585           82 :       CALL cp_fm_struct_release(fm_struct_tmp)
     586          246 :       ALLOCATE (preconditioner_env%full_evals(n))
     587          246 :       ALLOCATE (preconditioner_env%occ_evals(k))
     588              : 
     589              :       ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
     590              :       !    possibly shifted by an amount lambda,
     591              :       !    and the same spectrum as the original H matrix in the space orthogonal to the C0
     592              :       !    with P=C0 C0 ^ T
     593              :       !    (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
     594              :       !    we exploit that the C0 are already the ritz states of H
     595           82 :       CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
     596           82 :       CALL cp_fm_to_fm(matrix_c0, matrix_sc0)
     597           82 :       CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
     598           82 :       CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
     599              : 
     600              :       ! An aside, try to estimate the error on the ritz values, we'll need it later on
     601              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
     602              :                                context=preconditioner_env%ctxt, &
     603           82 :                                para_env=preconditioner_env%para_env)
     604           82 :       CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
     605           82 :       CALL cp_fm_struct_release(fm_struct_tmp)
     606              :       ! since we only use diagonal elements this is a bit of a waste
     607           82 :       CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1)
     608          164 :       ALLOCATE (diag(k))
     609           82 :       CALL cp_fm_get_diag(matrix_s1, diag)
     610          918 :       error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
     611           82 :       DEALLOCATE (diag)
     612           82 :       CALL cp_fm_release(matrix_s1)
     613              :       ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
     614              :       ! is small enough. A large error combined with a small energy gap would otherwise lead to
     615              :       ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
     616              :       ! aggressively
     617           82 :       preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
     618              : 
     619           82 :       matrix_pre => preconditioner_env%fm
     620           82 :       CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
     621           82 :       CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
     622              :       ! tmp = H ( 1 - PS )
     623           82 :       CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
     624              : 
     625              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
     626              :                                context=preconditioner_env%ctxt, &
     627           82 :                                para_env=preconditioner_env%para_env)
     628           82 :       CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
     629           82 :       CALL cp_fm_struct_release(fm_struct_tmp)
     630           82 :       CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
     631              :       ! tmp = (1 - PS)^T H (1-PS)
     632           82 :       CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
     633           82 :       CALL cp_fm_release(matrix_left)
     634              : 
     635          164 :       ALLOCATE (shifted_evals(k))
     636           82 :       lambda = lambda_base + error_estimate
     637          836 :       shifted_evals = c0_evals - lambda
     638           82 :       CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
     639           82 :       CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
     640           82 :       CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
     641              : 
     642              :       ! 2) diagonalize this operator
     643           82 :       CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
     644              : 
     645              :       ! test that the subspace remained conserved
     646              :       IF (.FALSE.) THEN
     647              :          CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
     648              :          CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
     649              :                                   context=preconditioner_env%ctxt, &
     650              :                                   para_env=preconditioner_env%para_env)
     651              :          CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
     652              :          CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
     653              :          CALL cp_fm_struct_release(fm_struct_tmp)
     654              :          ALLOCATE (norms(k))
     655              :          CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
     656              :          CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
     657              : 
     658              :          WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
     659              :          DEALLOCATE (norms)
     660              :          CALL cp_fm_release(matrix_s1)
     661              :          CALL cp_fm_release(matrix_s2)
     662              :       END IF
     663              : 
     664              :       ! 3) replace the lowest k evals and evecs with what they should be
     665          836 :       preconditioner_env%occ_evals = c0_evals
     666              :       ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all)
     667          836 :       preconditioner_env%full_evals(1:k) = c0_evals
     668           82 :       CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
     669              : 
     670           82 :       CALL cp_fm_release(matrix_sc0)
     671           82 :       CALL cp_fm_release(matrix_hc0)
     672           82 :       CALL cp_fm_release(matrix_tmp)
     673           82 :       DEALLOCATE (shifted_evals)
     674              : 
     675           82 :       CALL timestop(handle)
     676              : 
     677          574 :    END SUBROUTINE make_full_all_ortho
     678              : 
     679              : ! **************************************************************************************************
     680              : !> \brief generates a preconditioner matrix H-lambda S+(SC)(2.0*CT*H*C+delta)(SC)^T
     681              : !>        for later inversion.
     682              : !>        H is the Kohn Sham matrix
     683              : !>        lambda*S shifts the spectrum of the generalized form up by lambda
     684              : !>        the last term only shifts the occupied space (reversing them in energy order)
     685              : !>        This form is implicitly multiplied from both sides by S^0.5
     686              : !>        This ensures we precondition the correct quantity
     687              : !>        Before this reads S^-0.5 H S^-0.5 + lambda + (S^0.5 C)shifts(S^0.5 C)T
     688              : !>        which might be a bit more obvious
     689              : !>        Replaced the old full_single_inverse at revision 14616
     690              : !> \param preconditioner_env the preconditioner env
     691              : !> \param matrix_c0 the MO coefficient matrix (fm)
     692              : !> \param matrix_h Kohn-Sham matrix (dbcsr)
     693              : !> \param energy_gap an additional shift in lambda=-E_homo+energy_gap
     694              : !> \param matrix_s the overlap matrix if not orthonormal (dbcsr, optional)
     695              : ! **************************************************************************************************
     696         4045 :    SUBROUTINE make_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, energy_gap, matrix_s)
     697              :       TYPE(preconditioner_type)                          :: preconditioner_env
     698              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_c0
     699              :       TYPE(dbcsr_type), POINTER                          :: matrix_h
     700              :       REAL(KIND=dp)                                      :: energy_gap
     701              :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
     702              : 
     703              :       CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_inverse'
     704              : 
     705              :       INTEGER                                            :: handle, k, n
     706              :       REAL(KIND=dp)                                      :: max_ev, min_ev, pre_shift
     707              :       TYPE(arnoldi_env_type)                             :: arnoldi_env
     708         4045 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrices
     709              :       TYPE(dbcsr_type), TARGET                           :: dbcsr_cThc, dbcsr_hc, dbcsr_sc, mo_dbcsr
     710              : 
     711         4045 :       CALL timeset(routineN, handle)
     712              : 
     713              :       ! Allocate all working matrices needed
     714         4045 :       CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
     715              :       ! copy the fm MO's to a sparse matrix, can be solved better if the sparse version is already present
     716              :       ! but for the time beeing this will do
     717         4045 :       CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h)
     718         4045 :       CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr)
     719         4045 :       CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr)
     720         4045 :       CALL cp_dbcsr_m_by_n_from_template(dbcsr_cThc, matrix_h, k, k, sym=dbcsr_type_symmetric)
     721              : 
     722              :       ! Check whether the output matrix was already created, if not do it now
     723         4045 :       IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
     724         4045 :          ALLOCATE (preconditioner_env%sparse_matrix)
     725              :       END IF
     726              : 
     727              :       ! Put the first term of the preconditioner (H) into the output matrix
     728         4045 :       CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_h)
     729              : 
     730              :       ! Precompute some matrices
     731              :       ! S*C, if orthonormal this will be simply C so a copy will do
     732         4045 :       IF (PRESENT(matrix_s)) THEN
     733         3657 :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc)
     734              :       ELSE
     735          388 :          CALL dbcsr_copy(dbcsr_sc, mo_dbcsr)
     736              :       END IF
     737              : 
     738              : !----------------------------compute the occupied subspace and shift it ------------------------------------
     739              :       ! cT*H*C which will be used to shift the occupied states to 0
     740         4045 :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc)
     741         4045 :       CALL dbcsr_multiply("T", "N", 1.0_dp, mo_dbcsr, dbcsr_hc, 0.0_dp, dbcsr_cThc)
     742              : 
     743              :       ! Compute the Energy of the HOMO. We will use this as a reference energy
     744         8090 :       ALLOCATE (matrices(1))
     745         4045 :       matrices(1)%matrix => dbcsr_cThc
     746              :       CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=1.0E-3_dp, selection_crit=2, &
     747         4045 :                              nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
     748         4045 :       IF (ASSOCIATED(preconditioner_env%max_ev_vector)) &
     749         2310 :          CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%max_ev_vector)
     750         4045 :       CALL arnoldi_ev(matrices, arnoldi_env)
     751         4045 :       max_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
     752              : 
     753              :       ! save the ev as guess for the next time
     754         4045 :       IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector)
     755         4045 :       CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector)
     756         4045 :       CALL deallocate_arnoldi_env(arnoldi_env)
     757         4045 :       DEALLOCATE (matrices)
     758              : 
     759              :       ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H
     760         4045 :       CALL dbcsr_add_on_diag(dbcsr_cThc, -0.5_dp)
     761              :       ! Get the AO representation of the shift (see above why S is needed), W-matrix like object
     762         4045 :       CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cThc, 0.0_dp, dbcsr_hc)
     763         4045 :       CALL dbcsr_multiply("N", "T", -1.0_dp, dbcsr_hc, dbcsr_sc, 1.0_dp, preconditioner_env%sparse_matrix)
     764              : 
     765              : !-------------------------------------compute eigenvalues of H ----------------------------------------------
     766              :       ! Setup the arnoldi procedure to compute the lowest ev. if S is present this has to be the generalized ev
     767         4045 :       IF (PRESENT(matrix_s)) THEN
     768        10971 :          ALLOCATE (matrices(2))
     769         3657 :          matrices(1)%matrix => preconditioner_env%sparse_matrix
     770         3657 :          matrices(2)%matrix => matrix_s
     771              :          CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
     772         3657 :                                 nval_request=1, nrestarts=21, generalized_ev=.TRUE., iram=.FALSE.)
     773              :       ELSE
     774          776 :          ALLOCATE (matrices(1))
     775          388 :          matrices(1)%matrix => preconditioner_env%sparse_matrix
     776              :          CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
     777          388 :                                 nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
     778              :       END IF
     779         4045 :       IF (ASSOCIATED(preconditioner_env%min_ev_vector)) &
     780         2310 :          CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%min_ev_vector)
     781              : 
     782              :       ! compute the LUMO energy
     783         4045 :       CALL arnoldi_ev(matrices, arnoldi_env)
     784         4045 :       min_eV = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
     785              : 
     786              :       ! save the lumo vector for restarting in the next step
     787         4045 :       IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector)
     788         4045 :       CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector)
     789         4045 :       CALL deallocate_arnoldi_env(arnoldi_env)
     790         4045 :       DEALLOCATE (matrices)
     791              : 
     792              : !-------------------------------------compute eigenvalues of H ----------------------------------------------
     793              :       ! Shift the Lumo to the 1.5*the computed energy_gap or the external energy gap value
     794              :       ! The factor 1.5 is determined by trying. If the LUMO is positive, enough, just leave it alone
     795         4045 :       pre_shift = MAX(1.5_dp*(min_ev - max_ev), energy_gap)
     796         4045 :       IF (min_ev .LT. pre_shift) THEN
     797         4037 :          pre_shift = pre_shift - min_ev
     798              :       ELSE
     799            8 :          pre_shift = 0.0_dp
     800              :       END IF
     801         4045 :       IF (PRESENT(matrix_s)) THEN
     802         3657 :          CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift)
     803              :       ELSE
     804          388 :          CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift)
     805              :       END IF
     806              : 
     807         4045 :       CALL dbcsr_release(mo_dbcsr)
     808         4045 :       CALL dbcsr_release(dbcsr_hc)
     809         4045 :       CALL dbcsr_release(dbcsr_sc)
     810         4045 :       CALL dbcsr_release(dbcsr_cThc)
     811              : 
     812         4045 :       CALL timestop(handle)
     813              : 
     814         4045 :    END SUBROUTINE make_full_single_inverse
     815              : 
     816              : END MODULE preconditioner_makes
     817              : 
        

Generated by: LCOV version 2.0-1