LCOV - code coverage report
Current view: top level - src - qs_mo_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:c740a4b) Lines: 259 293 88.4 %
Date: 2025-05-29 08:02:39 Functions: 10 11 90.9 %

          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 collects routines that perform operations directly related to MOs
      10             : !> \note
      11             : !>      first version : most routines imported
      12             : !> \author Joost VandeVondele (2003-08)
      13             : ! **************************************************************************************************
      14             : MODULE qs_mo_methods
      15             :    USE admm_types,                      ONLY: admm_type
      16             :    USE admm_utils,                      ONLY: admm_correct_for_eigenvalues,&
      17             :                                               admm_uncorrect_for_eigenvalues
      18             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      19             :    USE cp_control_types,                ONLY: hairy_probes_type
      20             :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      21             :                                               dbcsr_get_info,&
      22             :                                               dbcsr_init_p,&
      23             :                                               dbcsr_multiply,&
      24             :                                               dbcsr_p_type,&
      25             :                                               dbcsr_release_p,&
      26             :                                               dbcsr_type,&
      27             :                                               dbcsr_type_no_symmetry
      28             :    USE cp_dbcsr_diag,                   ONLY: cp_dbcsr_syevd,&
      29             :                                               cp_dbcsr_syevx
      30             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      31             :                                               copy_fm_to_dbcsr,&
      32             :                                               cp_dbcsr_m_by_n_from_template,&
      33             :                                               cp_dbcsr_sm_fm_multiply
      34             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_syrk,&
      35             :                                               cp_fm_triangular_multiply
      36             :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose
      37             :    USE cp_fm_diag,                      ONLY: choose_eigv_solver,&
      38             :                                               cp_fm_power,&
      39             :                                               cp_fm_syevx
      40             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      41             :                                               cp_fm_struct_release,&
      42             :                                               cp_fm_struct_type
      43             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      44             :                                               cp_fm_get_info,&
      45             :                                               cp_fm_release,&
      46             :                                               cp_fm_to_fm,&
      47             :                                               cp_fm_type
      48             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      49             :    USE kinds,                           ONLY: dp
      50             :    USE message_passing,                 ONLY: mp_para_env_type
      51             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      52             :    USE physcon,                         ONLY: evolt
      53             :    USE qs_mo_occupation,                ONLY: set_mo_occupation
      54             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      55             :                                               mo_set_type
      56             :    USE scf_control_types,               ONLY: scf_control_type
      57             : #include "./base/base_uses.f90"
      58             : 
      59             :    IMPLICIT NONE
      60             :    PRIVATE
      61             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_methods'
      62             : 
      63             :    PUBLIC :: make_basis_simple, make_basis_cholesky, make_basis_sv, make_basis_sm, &
      64             :              make_basis_lowdin, calculate_subspace_eigenvalues, &
      65             :              calculate_orthonormality, calculate_magnitude, make_mo_eig
      66             : 
      67             :    INTERFACE calculate_subspace_eigenvalues
      68             :       MODULE PROCEDURE subspace_eigenvalues_ks_fm
      69             :       MODULE PROCEDURE subspace_eigenvalues_ks_dbcsr
      70             :    END INTERFACE
      71             : 
      72             :    INTERFACE make_basis_sv
      73             :       MODULE PROCEDURE make_basis_sv_fm
      74             :       MODULE PROCEDURE make_basis_sv_dbcsr
      75             :    END INTERFACE
      76             : 
      77             : CONTAINS
      78             : 
      79             : ! **************************************************************************************************
      80             : !> \brief returns an S-orthonormal basis v (v^T S v ==1)
      81             : !> \param vmatrix ...
      82             : !> \param ncol ...
      83             : !> \param matrix_s ...
      84             : !> \par History
      85             : !>      03.2006 created [Joost VandeVondele]
      86             : ! **************************************************************************************************
      87       33118 :    SUBROUTINE make_basis_sm(vmatrix, ncol, matrix_s)
      88             :       TYPE(cp_fm_type), INTENT(IN)                       :: vmatrix
      89             :       INTEGER, INTENT(IN)                                :: ncol
      90             :       TYPE(dbcsr_type), POINTER                          :: matrix_s
      91             : 
      92             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'make_basis_sm'
      93             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
      94             : 
      95             :       INTEGER                                            :: handle, n, ncol_global
      96             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
      97             :       TYPE(cp_fm_type)                                   :: overlap_vv, svmatrix
      98             : 
      99         128 :       IF (ncol .EQ. 0) RETURN
     100             : 
     101        6598 :       CALL timeset(routineN, handle)
     102             : 
     103        6598 :       CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
     104        6598 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     105             : 
     106        6598 :       CALL cp_fm_create(svmatrix, vmatrix%matrix_struct, "SV")
     107        6598 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, svmatrix, ncol)
     108             : 
     109        6598 :       NULLIFY (fm_struct_tmp)
     110             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     111             :                                para_env=vmatrix%matrix_struct%para_env, &
     112        6598 :                                context=vmatrix%matrix_struct%context)
     113        6598 :       CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
     114        6598 :       CALL cp_fm_struct_release(fm_struct_tmp)
     115             : 
     116        6598 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
     117        6598 :       CALL cp_fm_cholesky_decompose(overlap_vv)
     118        6598 :       CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
     119             : 
     120        6598 :       CALL cp_fm_release(overlap_vv)
     121        6598 :       CALL cp_fm_release(svmatrix)
     122             : 
     123        6598 :       CALL timestop(handle)
     124             : 
     125        6726 :    END SUBROUTINE make_basis_sm
     126             : 
     127             : ! **************************************************************************************************
     128             : !> \brief returns an S-orthonormal basis v and the corresponding matrix S*v as well
     129             : !> \param vmatrix ...
     130             : !> \param ncol ...
     131             : !> \param svmatrix ...
     132             : !> \par History
     133             : !>      03.2006 created [Joost VandeVondele]
     134             : ! **************************************************************************************************
     135           0 :    SUBROUTINE make_basis_sv_fm(vmatrix, ncol, svmatrix)
     136             : 
     137             :       TYPE(cp_fm_type), INTENT(IN)                       :: vmatrix
     138             :       INTEGER, INTENT(IN)                                :: ncol
     139             :       TYPE(cp_fm_type), INTENT(IN)                       :: svmatrix
     140             : 
     141             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'make_basis_sv_fm'
     142             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     143             : 
     144             :       INTEGER                                            :: handle, n, ncol_global
     145             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     146             :       TYPE(cp_fm_type)                                   :: overlap_vv
     147             : 
     148           0 :       IF (ncol .EQ. 0) RETURN
     149             : 
     150           0 :       CALL timeset(routineN, handle)
     151           0 :       NULLIFY (fm_struct_tmp)
     152             : 
     153           0 :       CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
     154           0 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     155             : 
     156             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     157             :                                para_env=vmatrix%matrix_struct%para_env, &
     158           0 :                                context=vmatrix%matrix_struct%context)
     159           0 :       CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
     160           0 :       CALL cp_fm_struct_release(fm_struct_tmp)
     161             : 
     162           0 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
     163           0 :       CALL cp_fm_cholesky_decompose(overlap_vv)
     164           0 :       CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
     165           0 :       CALL cp_fm_triangular_multiply(overlap_vv, svmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
     166             : 
     167           0 :       CALL cp_fm_release(overlap_vv)
     168             : 
     169           0 :       CALL timestop(handle)
     170             : 
     171           0 :    END SUBROUTINE make_basis_sv_fm
     172             : 
     173             : ! **************************************************************************************************
     174             : !> \brief ...
     175             : !> \param vmatrix ...
     176             : !> \param ncol ...
     177             : !> \param svmatrix ...
     178             : !> \param para_env ...
     179             : !> \param blacs_env ...
     180             : ! **************************************************************************************************
     181        2834 :    SUBROUTINE make_basis_sv_dbcsr(vmatrix, ncol, svmatrix, para_env, blacs_env)
     182             : 
     183             :       TYPE(dbcsr_type)                                   :: vmatrix
     184             :       INTEGER, INTENT(IN)                                :: ncol
     185             :       TYPE(dbcsr_type)                                   :: svmatrix
     186             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     187             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     188             : 
     189             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_sv_dbcsr'
     190             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     191             : 
     192             :       INTEGER                                            :: handle, n, ncol_global
     193             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     194             :       TYPE(cp_fm_type)                                   :: fm_svmatrix, fm_vmatrix, overlap_vv
     195             : 
     196          74 :       IF (ncol .EQ. 0) RETURN
     197             : 
     198         460 :       CALL timeset(routineN, handle)
     199             : 
     200             :       !CALL cp_fm_get_info(matrix=vmatrix,nrow_global=n,ncol_global=ncol_global)
     201         460 :       CALL dbcsr_get_info(vmatrix, nfullrows_total=n, nfullcols_total=ncol_global)
     202         460 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     203             : 
     204             :       CALL cp_fm_struct_create(fm_struct_tmp, context=blacs_env, nrow_global=ncol, &
     205         460 :                                ncol_global=ncol, para_env=para_env)
     206         460 :       CALL cp_fm_create(overlap_vv, fm_struct_tmp, name="fm_overlap_vv")
     207         460 :       CALL cp_fm_struct_release(fm_struct_tmp)
     208             : 
     209             :       CALL cp_fm_struct_create(fm_struct_tmp, context=blacs_env, nrow_global=n, &
     210         460 :                                ncol_global=ncol_global, para_env=para_env)
     211         460 :       CALL cp_fm_create(fm_vmatrix, fm_struct_tmp, name="fm_vmatrix")
     212         460 :       CALL cp_fm_create(fm_svmatrix, fm_struct_tmp, name="fm_svmatrix")
     213         460 :       CALL cp_fm_struct_release(fm_struct_tmp)
     214             : 
     215         460 :       CALL copy_dbcsr_to_fm(vmatrix, fm_vmatrix)
     216         460 :       CALL copy_dbcsr_to_fm(svmatrix, fm_svmatrix)
     217             : 
     218         460 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, fm_vmatrix, fm_svmatrix, rzero, overlap_vv)
     219         460 :       CALL cp_fm_cholesky_decompose(overlap_vv)
     220         460 :       CALL cp_fm_triangular_multiply(overlap_vv, fm_vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
     221         460 :       CALL cp_fm_triangular_multiply(overlap_vv, fm_svmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
     222             : 
     223         460 :       CALL copy_fm_to_dbcsr(fm_vmatrix, vmatrix)
     224         460 :       CALL copy_fm_to_dbcsr(fm_svmatrix, svmatrix)
     225             : 
     226         460 :       CALL cp_fm_release(overlap_vv)
     227         460 :       CALL cp_fm_release(fm_vmatrix)
     228         460 :       CALL cp_fm_release(fm_svmatrix)
     229             : 
     230         460 :       CALL timestop(handle)
     231             : 
     232         534 :    END SUBROUTINE make_basis_sv_dbcsr
     233             : 
     234             : ! **************************************************************************************************
     235             : !> \brief return a set of S orthonormal vectors (C^T S C == 1) where
     236             : !>      the cholesky decomposed form of S is passed as an argument
     237             : !> \param vmatrix ...
     238             : !> \param ncol ...
     239             : !> \param ortho cholesky decomposed S matrix
     240             : !> \par History
     241             : !>      03.2006 created [Joost VandeVondele]
     242             : !> \note
     243             : !>      if the cholesky decomposed S matrix is not available
     244             : !>      use make_basis_sm since this is much faster than computing the
     245             : !>      cholesky decomposition of S
     246             : ! **************************************************************************************************
     247       24416 :    SUBROUTINE make_basis_cholesky(vmatrix, ncol, ortho)
     248             : 
     249             :       TYPE(cp_fm_type), INTENT(IN)                       :: vmatrix
     250             :       INTEGER, INTENT(IN)                                :: ncol
     251             :       TYPE(cp_fm_type), INTENT(IN)                       :: ortho
     252             : 
     253             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_cholesky'
     254             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     255             : 
     256             :       INTEGER                                            :: handle, n, ncol_global
     257             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     258             :       TYPE(cp_fm_type)                                   :: overlap_vv
     259             : 
     260          80 :       IF (ncol .EQ. 0) RETURN
     261             : 
     262        6084 :       CALL timeset(routineN, handle)
     263        6084 :       NULLIFY (fm_struct_tmp)
     264             : 
     265        6084 :       CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
     266        6084 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     267             : 
     268             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     269             :                                para_env=vmatrix%matrix_struct%para_env, &
     270        6084 :                                context=vmatrix%matrix_struct%context)
     271        6084 :       CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
     272        6084 :       CALL cp_fm_struct_release(fm_struct_tmp)
     273             : 
     274        6084 :       CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol)
     275        6084 :       CALL cp_fm_syrk('U', 'T', n, rone, vmatrix, 1, 1, rzero, overlap_vv)
     276        6084 :       CALL cp_fm_cholesky_decompose(overlap_vv)
     277        6084 :       CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
     278        6084 :       CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol, invert_tr=.TRUE.)
     279             : 
     280        6084 :       CALL cp_fm_release(overlap_vv)
     281             : 
     282        6084 :       CALL timestop(handle)
     283             : 
     284        6164 :    END SUBROUTINE make_basis_cholesky
     285             : 
     286             : ! **************************************************************************************************
     287             : !> \brief return a set of S orthonormal vectors (C^T S C == 1) where
     288             : !>      a Loedwin transformation is applied to keep the rotated vectors as close
     289             : !>      as possible to the original ones
     290             : !> \param vmatrix ...
     291             : !> \param ncol ...
     292             : !> \param matrix_s ...
     293             : !> \param
     294             : !> \par History
     295             : !>      05.2009 created [MI]
     296             : !> \note
     297             : ! **************************************************************************************************
     298        2772 :    SUBROUTINE make_basis_lowdin(vmatrix, ncol, matrix_s)
     299             : 
     300             :       TYPE(cp_fm_type), INTENT(IN)                       :: vmatrix
     301             :       INTEGER, INTENT(IN)                                :: ncol
     302             :       TYPE(dbcsr_type)                                   :: matrix_s
     303             : 
     304             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'make_basis_lowdin'
     305             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     306             : 
     307             :       INTEGER                                            :: handle, n, ncol_global, ndep
     308             :       REAL(dp)                                           :: threshold
     309             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     310             :       TYPE(cp_fm_type)                                   :: csc, sc, work
     311             : 
     312           0 :       IF (ncol .EQ. 0) RETURN
     313             : 
     314         462 :       CALL timeset(routineN, handle)
     315         462 :       NULLIFY (fm_struct_tmp)
     316         462 :       threshold = 1.0E-7_dp
     317         462 :       CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
     318         462 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     319             : 
     320         462 :       CALL cp_fm_create(sc, vmatrix%matrix_struct, "SC")
     321         462 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
     322             : 
     323         462 :       NULLIFY (fm_struct_tmp)
     324             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     325             :                                para_env=vmatrix%matrix_struct%para_env, &
     326         462 :                                context=vmatrix%matrix_struct%context)
     327         462 :       CALL cp_fm_create(csc, fm_struct_tmp, "csc")
     328         462 :       CALL cp_fm_create(work, fm_struct_tmp, "work")
     329         462 :       CALL cp_fm_struct_release(fm_struct_tmp)
     330             : 
     331         462 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
     332         462 :       CALL cp_fm_power(csc, work, -0.5_dp, threshold, ndep)
     333         462 :       CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
     334         462 :       CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
     335             : 
     336         462 :       CALL cp_fm_release(csc)
     337         462 :       CALL cp_fm_release(sc)
     338         462 :       CALL cp_fm_release(work)
     339             : 
     340         462 :       CALL timestop(handle)
     341             : 
     342         462 :    END SUBROUTINE make_basis_lowdin
     343             : 
     344             : ! **************************************************************************************************
     345             : !> \brief given a set of vectors, return an orthogonal (C^T C == 1) set
     346             : !>      spanning the same space (notice, only for cases where S==1)
     347             : !> \param vmatrix ...
     348             : !> \param ncol ...
     349             : !> \par History
     350             : !>      03.2006 created [Joost VandeVondele]
     351             : ! **************************************************************************************************
     352       14526 :    SUBROUTINE make_basis_simple(vmatrix, ncol)
     353             : 
     354             :       TYPE(cp_fm_type), INTENT(IN)                       :: vmatrix
     355             :       INTEGER, INTENT(IN)                                :: ncol
     356             : 
     357             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'make_basis_simple'
     358             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     359             : 
     360             :       INTEGER                                            :: handle, n, ncol_global
     361             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     362             :       TYPE(cp_fm_type)                                   :: overlap_vv
     363             : 
     364          70 :       IF (ncol .EQ. 0) RETURN
     365             : 
     366        3614 :       CALL timeset(routineN, handle)
     367             : 
     368        3614 :       NULLIFY (fm_struct_tmp)
     369             : 
     370        3614 :       CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
     371        3614 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     372             : 
     373             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     374             :                                para_env=vmatrix%matrix_struct%para_env, &
     375        3614 :                                context=vmatrix%matrix_struct%context)
     376        3614 :       CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
     377        3614 :       CALL cp_fm_struct_release(fm_struct_tmp)
     378             : 
     379        3614 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, vmatrix, rzero, overlap_vv)
     380        3614 :       CALL cp_fm_cholesky_decompose(overlap_vv)
     381        3614 :       CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
     382             : 
     383        3614 :       CALL cp_fm_release(overlap_vv)
     384             : 
     385        3614 :       CALL timestop(handle)
     386             : 
     387        3684 :    END SUBROUTINE make_basis_simple
     388             : 
     389             : ! **************************************************************************************************
     390             : !> \brief computes ritz values of a set of orbitals given a ks_matrix
     391             : !>      rotates the orbitals into eigenstates depending on do_rotation
     392             : !>      writes the evals to the screen depending on ionode/scr
     393             : !> \param orbitals S-orthonormal orbitals
     394             : !> \param ks_matrix Kohn-Sham matrix
     395             : !> \param evals_arg optional, filled with the evals
     396             : !> \param ionode , scr if present write to unit scr where ionode
     397             : !> \param scr ...
     398             : !> \param do_rotation optional rotate orbitals (default=.TRUE.)
     399             : !>        note that rotating the orbitals is slower
     400             : !> \param co_rotate an optional set of orbitals rotated by the same rotation matrix
     401             : !> \param co_rotate_dbcsr ...
     402             : !> \par History
     403             : !>      08.2004 documented and added do_rotation [Joost VandeVondele]
     404             : !>      09.2008 only compute eigenvalues if rotation is not needed
     405             : ! **************************************************************************************************
     406         886 :    SUBROUTINE subspace_eigenvalues_ks_fm(orbitals, ks_matrix, evals_arg, ionode, scr, &
     407             :                                          do_rotation, co_rotate, co_rotate_dbcsr)
     408             : 
     409             :       TYPE(cp_fm_type), INTENT(IN)                       :: orbitals
     410             :       TYPE(dbcsr_type), POINTER                          :: ks_matrix
     411             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: evals_arg
     412             :       LOGICAL, INTENT(IN), OPTIONAL                      :: ionode
     413             :       INTEGER, INTENT(IN), OPTIONAL                      :: scr
     414             :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_rotation
     415             :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: co_rotate
     416             :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: co_rotate_dbcsr
     417             : 
     418             :       CHARACTER(len=*), PARAMETER :: routineN = 'subspace_eigenvalues_ks_fm'
     419             : 
     420             :       INTEGER                                            :: handle, i, j, n, ncol_global, nrow_global
     421             :       LOGICAL                                            :: compute_evecs, do_rotation_local
     422         886 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
     423             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     424             :       TYPE(cp_fm_type)                                   :: e_vectors, h_block, weighted_vectors, &
     425             :                                                             weighted_vectors2
     426             : 
     427         886 :       CALL timeset(routineN, handle)
     428             : 
     429         886 :       do_rotation_local = .TRUE.
     430         886 :       IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
     431             : 
     432         886 :       NULLIFY (fm_struct_tmp)
     433             :       CALL cp_fm_get_info(matrix=orbitals, &
     434             :                           ncol_global=ncol_global, &
     435         886 :                           nrow_global=nrow_global)
     436             : 
     437             :       IF (do_rotation_local) THEN
     438             :          compute_evecs = .TRUE.
     439             :       ELSE
     440             :          ! this would be the logical choice if syevx computing only evals were faster than syevd computing evecs and evals.
     441             :          compute_evecs = .FALSE.
     442             :          ! this is not the case, so lets compute evecs always
     443             :          compute_evecs = .TRUE.
     444             :       END IF
     445             : 
     446         886 :       IF (ncol_global .GT. 0) THEN
     447             : 
     448        2436 :          ALLOCATE (evals(ncol_global))
     449             : 
     450         812 :          CALL cp_fm_create(weighted_vectors, orbitals%matrix_struct, "weighted_vectors")
     451             :          CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, &
     452             :                                   para_env=orbitals%matrix_struct%para_env, &
     453         812 :                                   context=orbitals%matrix_struct%context)
     454         812 :          CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
     455         812 :          IF (compute_evecs) THEN
     456         812 :             CALL cp_fm_create(e_vectors, fm_struct_tmp, name="e vectors")
     457             :          END IF
     458         812 :          CALL cp_fm_struct_release(fm_struct_tmp)
     459             : 
     460             :          ! h subblock and diag
     461         812 :          CALL cp_dbcsr_sm_fm_multiply(ks_matrix, orbitals, weighted_vectors, ncol_global)
     462             : 
     463             :          CALL parallel_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 1.0_dp, &
     464         812 :                             orbitals, weighted_vectors, 0.0_dp, h_block)
     465             : 
     466             :          ! if eigenvectors are required, go for syevd, otherwise only compute eigenvalues
     467             :          IF (compute_evecs) THEN
     468         812 :             CALL choose_eigv_solver(h_block, e_vectors, evals)
     469             :          ELSE
     470             :             CALL cp_fm_syevx(h_block, eigenvalues=evals)
     471             :          END IF
     472             : 
     473             :          ! rotate the orbitals
     474         812 :          IF (do_rotation_local) THEN
     475             :             CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
     476         528 :                                orbitals, e_vectors, 0.0_dp, weighted_vectors)
     477         528 :             CALL cp_fm_to_fm(weighted_vectors, orbitals)
     478         528 :             IF (PRESENT(co_rotate)) THEN
     479             :                CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
     480           0 :                                   co_rotate, e_vectors, 0.0_dp, weighted_vectors)
     481           0 :                CALL cp_fm_to_fm(weighted_vectors, co_rotate)
     482             :             END IF
     483         528 :             IF (PRESENT(co_rotate_dbcsr)) THEN
     484         136 :                IF (ASSOCIATED(co_rotate_dbcsr)) THEN
     485         136 :                   CALL cp_fm_create(weighted_vectors2, orbitals%matrix_struct, "weighted_vectors")
     486         136 :                   CALL copy_dbcsr_to_fm(co_rotate_dbcsr, weighted_vectors2)
     487             :                   CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
     488         136 :                                      weighted_vectors2, e_vectors, 0.0_dp, weighted_vectors)
     489         136 :                   CALL copy_fm_to_dbcsr(weighted_vectors, co_rotate_dbcsr)
     490         136 :                   CALL cp_fm_release(weighted_vectors2)
     491             :                END IF
     492             :             END IF
     493             :          END IF
     494             : 
     495             :          ! give output
     496         812 :          IF (PRESENT(evals_arg)) THEN
     497         780 :             n = MIN(SIZE(evals_arg), SIZE(evals))
     498        7202 :             evals_arg(1:n) = evals(1:n)
     499             :          END IF
     500             : 
     501         812 :          IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
     502         180 :             IF (.NOT. PRESENT(ionode)) CPABORT("IONODE?")
     503         180 :             IF (.NOT. PRESENT(scr)) CPABORT("SCR?")
     504         180 :             IF (ionode) THEN
     505         225 :                DO i = 1, ncol_global, 4
     506         135 :                   j = MIN(3, ncol_global - i)
     507          90 :                   SELECT CASE (j)
     508             :                   CASE (3)
     509         104 :                      WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
     510             :                   CASE (2)
     511           2 :                      WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
     512             :                   CASE (1)
     513           7 :                      WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
     514             :                   CASE (0)
     515         135 :                      WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
     516             :                   END SELECT
     517             :                END DO
     518             :             END IF
     519             :          END IF
     520             : 
     521         812 :          CALL cp_fm_release(weighted_vectors)
     522         812 :          CALL cp_fm_release(h_block)
     523             :          IF (compute_evecs) THEN
     524         812 :             CALL cp_fm_release(e_vectors)
     525             :          END IF
     526             : 
     527        2436 :          DEALLOCATE (evals)
     528             : 
     529             :       END IF
     530             : 
     531         886 :       CALL timestop(handle)
     532             : 
     533        1772 :    END SUBROUTINE subspace_eigenvalues_ks_fm
     534             : 
     535             : ! **************************************************************************************************
     536             : !> \brief ...
     537             : !> \param orbitals ...
     538             : !> \param ks_matrix ...
     539             : !> \param evals_arg ...
     540             : !> \param ionode ...
     541             : !> \param scr ...
     542             : !> \param do_rotation ...
     543             : !> \param co_rotate ...
     544             : !> \param para_env ...
     545             : !> \param blacs_env ...
     546             : ! **************************************************************************************************
     547        5328 :    SUBROUTINE subspace_eigenvalues_ks_dbcsr(orbitals, ks_matrix, evals_arg, ionode, scr, &
     548             :                                             do_rotation, co_rotate, para_env, blacs_env)
     549             : 
     550             :       TYPE(dbcsr_type), POINTER                          :: orbitals, ks_matrix
     551             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: evals_arg
     552             :       LOGICAL, INTENT(IN), OPTIONAL                      :: ionode
     553             :       INTEGER, INTENT(IN), OPTIONAL                      :: scr
     554             :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_rotation
     555             :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: co_rotate
     556             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     557             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     558             : 
     559             :       CHARACTER(len=*), PARAMETER :: routineN = 'subspace_eigenvalues_ks_dbcsr'
     560             : 
     561             :       INTEGER                                            :: handle, i, j, ncol_global, nrow_global
     562             :       LOGICAL                                            :: compute_evecs, do_rotation_local
     563        5328 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
     564             :       TYPE(dbcsr_type), POINTER                          :: e_vectors, h_block, weighted_vectors
     565             : 
     566        5328 :       CALL timeset(routineN, handle)
     567             : 
     568        5328 :       do_rotation_local = .TRUE.
     569        5328 :       IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
     570             : 
     571        5328 :       NULLIFY (e_vectors, h_block, weighted_vectors)
     572             : 
     573             :       CALL dbcsr_get_info(matrix=orbitals, &
     574             :                           nfullcols_total=ncol_global, &
     575        5328 :                           nfullrows_total=nrow_global)
     576             : 
     577             :       IF (do_rotation_local) THEN
     578             :          compute_evecs = .TRUE.
     579             :       ELSE
     580             :          ! this would be the logical choice if syevx computing only evals were faster than syevd computing evecs and evals.
     581             :          compute_evecs = .FALSE.
     582             :          ! this is not the case, so lets compute evecs always
     583             :          compute_evecs = .TRUE.
     584             :       END IF
     585             : 
     586        5328 :       IF (ncol_global .GT. 0) THEN
     587             : 
     588       15348 :          ALLOCATE (evals(ncol_global))
     589             : 
     590        5116 :          CALL dbcsr_init_p(weighted_vectors)
     591        5116 :          CALL dbcsr_copy(weighted_vectors, orbitals, name="weighted_vectors")
     592             : 
     593        5116 :          CALL dbcsr_init_p(h_block)
     594             :          CALL cp_dbcsr_m_by_n_from_template(h_block, template=orbitals, m=ncol_global, n=ncol_global, &
     595        5116 :                                             sym=dbcsr_type_no_symmetry)
     596             : 
     597             : !!!!!!!!!!!!!!XXXXXXXXXXXXXXXXXXX!!!!!!!!!!!!!
     598             :          IF (compute_evecs) THEN
     599        5116 :             CALL dbcsr_init_p(e_vectors)
     600             :             CALL cp_dbcsr_m_by_n_from_template(e_vectors, template=orbitals, m=ncol_global, n=ncol_global, &
     601        5116 :                                                sym=dbcsr_type_no_symmetry)
     602             :          END IF
     603             : 
     604             :          ! h subblock and diag
     605             :          CALL dbcsr_multiply('N', 'N', 1.0_dp, ks_matrix, orbitals, &
     606        5116 :                              0.0_dp, weighted_vectors)
     607             :          !CALL cp_dbcsr_sm_fm_multiply(ks_matrix,orbitals,weighted_vectors, ncol_global)
     608             : 
     609        5116 :          CALL dbcsr_multiply('T', 'N', 1.0_dp, orbitals, weighted_vectors, 0.0_dp, h_block)
     610             :          !CALL parallel_gemm('T','N',ncol_global,ncol_global,nrow_global,1.0_dp, &
     611             :          !                orbitals,weighted_vectors,0.0_dp,h_block)
     612             : 
     613             :          ! if eigenvectors are required, go for syevd, otherwise only compute eigenvalues
     614             :          IF (compute_evecs) THEN
     615        5116 :             CALL cp_dbcsr_syevd(h_block, e_vectors, evals, para_env=para_env, blacs_env=blacs_env)
     616             :          ELSE
     617             :             CALL cp_dbcsr_syevx(h_block, eigenvalues=evals, para_env=para_env, blacs_env=blacs_env)
     618             :          END IF
     619             : 
     620             :          ! rotate the orbitals
     621        5116 :          IF (do_rotation_local) THEN
     622        2504 :             CALL dbcsr_multiply('N', 'N', 1.0_dp, orbitals, e_vectors, 0.0_dp, weighted_vectors)
     623             :             !CALL parallel_gemm('N','N',nrow_global,ncol_global,ncol_global,1.0_dp, &
     624             :             !             orbitals,e_vectors,0.0_dp,weighted_vectors)
     625        2504 :             CALL dbcsr_copy(orbitals, weighted_vectors)
     626             :             !CALL cp_fm_to_fm(weighted_vectors,orbitals)
     627        2504 :             IF (PRESENT(co_rotate)) THEN
     628        2486 :                IF (ASSOCIATED(co_rotate)) THEN
     629        2486 :                   CALL dbcsr_multiply('N', 'N', 1.0_dp, co_rotate, e_vectors, 0.0_dp, weighted_vectors)
     630             :                   !CALL parallel_gemm('N','N',nrow_global,ncol_global,ncol_global,1.0_dp, &
     631             :                   !     co_rotate,e_vectors,0.0_dp,weighted_vectors)
     632        2486 :                   CALL dbcsr_copy(co_rotate, weighted_vectors)
     633             :                   !CALL cp_fm_to_fm(weighted_vectors,co_rotate)
     634             :                END IF
     635             :             END IF
     636             :          END IF
     637             : 
     638             :          ! give output
     639        5116 :          IF (PRESENT(evals_arg)) THEN
     640       14144 :             evals_arg(:) = evals(:)
     641             :          END IF
     642             : 
     643        5116 :          IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
     644           0 :             IF (.NOT. PRESENT(ionode)) CPABORT("IONODE?")
     645           0 :             IF (.NOT. PRESENT(scr)) CPABORT("SCR?")
     646           0 :             IF (ionode) THEN
     647           0 :                DO i = 1, ncol_global, 4
     648           0 :                   j = MIN(3, ncol_global - i)
     649           0 :                   SELECT CASE (j)
     650             :                   CASE (3)
     651           0 :                      WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
     652             :                   CASE (2)
     653           0 :                      WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
     654             :                   CASE (1)
     655           0 :                      WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
     656             :                   CASE (0)
     657           0 :                      WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
     658             :                   END SELECT
     659             :                END DO
     660             :             END IF
     661             :          END IF
     662             : 
     663        5116 :          CALL dbcsr_release_p(weighted_vectors)
     664        5116 :          CALL dbcsr_release_p(h_block)
     665             :          IF (compute_evecs) THEN
     666        5116 :             CALL dbcsr_release_p(e_vectors)
     667             :          END IF
     668             : 
     669        5116 :          DEALLOCATE (evals)
     670             : 
     671             :       END IF
     672             : 
     673        5328 :       CALL timestop(handle)
     674             : 
     675        5328 :    END SUBROUTINE subspace_eigenvalues_ks_dbcsr
     676             : 
     677             : ! computes the effective orthonormality of a set of mos given an s-matrix
     678             : ! orthonormality is the max deviation from unity of the C^T S C
     679             : ! **************************************************************************************************
     680             : !> \brief ...
     681             : !> \param orthonormality ...
     682             : !> \param mo_array ...
     683             : !> \param matrix_s ...
     684             : ! **************************************************************************************************
     685        1332 :    SUBROUTINE calculate_orthonormality(orthonormality, mo_array, matrix_s)
     686             :       REAL(KIND=dp)                                      :: orthonormality
     687             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     688             :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
     689             : 
     690             :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_orthonormality'
     691             : 
     692             :       INTEGER                                            :: handle, i, ispin, j, k, n, ncol_local, &
     693             :                                                             nrow_local, nspin
     694        1332 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     695             :       REAL(KIND=dp)                                      :: alpha, max_alpha
     696             :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
     697             :       TYPE(cp_fm_type)                                   :: overlap, svec
     698             : 
     699        1332 :       NULLIFY (tmp_fm_struct)
     700             : 
     701        1332 :       CALL timeset(routineN, handle)
     702             : 
     703        1332 :       nspin = SIZE(mo_array)
     704        1332 :       max_alpha = 0.0_dp
     705             : 
     706        2676 :       DO ispin = 1, nspin
     707        1344 :          IF (PRESENT(matrix_s)) THEN
     708             :             ! get S*C
     709        1286 :             CALL cp_fm_create(svec, mo_array(ispin)%mo_coeff%matrix_struct)
     710             :             CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
     711        1286 :                                 nrow_global=n, ncol_global=k)
     712             :             CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_array(ispin)%mo_coeff, &
     713        1286 :                                          svec, k)
     714             :             ! get C^T (S*C)
     715             :             CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
     716             :                                      para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
     717        1286 :                                      context=mo_array(ispin)%mo_coeff%matrix_struct%context)
     718        1286 :             CALL cp_fm_create(overlap, tmp_fm_struct)
     719        1286 :             CALL cp_fm_struct_release(tmp_fm_struct)
     720             :             CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
     721        1286 :                                svec, 0.0_dp, overlap)
     722        1286 :             CALL cp_fm_release(svec)
     723             :          ELSE
     724             :             ! orthogonal basis C^T C
     725             :             CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
     726          58 :                                 nrow_global=n, ncol_global=k)
     727             :             CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
     728             :                                      para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
     729          58 :                                      context=mo_array(ispin)%mo_coeff%matrix_struct%context)
     730          58 :             CALL cp_fm_create(overlap, tmp_fm_struct)
     731          58 :             CALL cp_fm_struct_release(tmp_fm_struct)
     732             :             CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
     733          58 :                                mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
     734             :          END IF
     735             :          CALL cp_fm_get_info(overlap, nrow_local=nrow_local, ncol_local=ncol_local, &
     736        1344 :                              row_indices=row_indices, col_indices=col_indices)
     737        7678 :          DO i = 1, nrow_local
     738      359494 :             DO j = 1, ncol_local
     739      351816 :                alpha = overlap%local_data(i, j)
     740      351816 :                IF (row_indices(i) .EQ. col_indices(j)) alpha = alpha - 1.0_dp
     741      358150 :                max_alpha = MAX(max_alpha, ABS(alpha))
     742             :             END DO
     743             :          END DO
     744        4020 :          CALL cp_fm_release(overlap)
     745             :       END DO
     746        1332 :       CALL mo_array(1)%mo_coeff%matrix_struct%para_env%max(max_alpha)
     747        1332 :       orthonormality = max_alpha
     748             :       ! write(*,*) "max deviation from orthonormalization ",orthonormality
     749             : 
     750        1332 :       CALL timestop(handle)
     751             : 
     752        1332 :    END SUBROUTINE calculate_orthonormality
     753             : 
     754             : ! computes the minimum/maximum magnitudes of C^T C. This could be useful
     755             : ! to detect problems in the case of nearly singular overlap matrices.
     756             : ! in this case, we expect the ratio of min/max to be large
     757             : ! this routine is only similar to mo_orthonormality if S==1
     758             : ! **************************************************************************************************
     759             : !> \brief ...
     760             : !> \param mo_array ...
     761             : !> \param mo_mag_min ...
     762             : !> \param mo_mag_max ...
     763             : ! **************************************************************************************************
     764        1312 :    SUBROUTINE calculate_magnitude(mo_array, mo_mag_min, mo_mag_max)
     765             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     766             :       REAL(KIND=dp)                                      :: mo_mag_min, mo_mag_max
     767             : 
     768             :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_magnitude'
     769             : 
     770             :       INTEGER                                            :: handle, ispin, k, n, nspin
     771        1312 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
     772             :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
     773             :       TYPE(cp_fm_type)                                   :: evecs, overlap
     774             : 
     775        1312 :       NULLIFY (tmp_fm_struct)
     776             : 
     777        1312 :       CALL timeset(routineN, handle)
     778             : 
     779        1312 :       nspin = SIZE(mo_array)
     780        1312 :       mo_mag_min = HUGE(0.0_dp)
     781        1312 :       mo_mag_max = -HUGE(0.0_dp)
     782        2624 :       DO ispin = 1, nspin
     783             :          CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
     784        1312 :                              nrow_global=n, ncol_global=k)
     785        3936 :          ALLOCATE (evals(k))
     786             :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
     787             :                                   para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
     788        1312 :                                   context=mo_array(ispin)%mo_coeff%matrix_struct%context)
     789        1312 :          CALL cp_fm_create(overlap, tmp_fm_struct)
     790        1312 :          CALL cp_fm_create(evecs, tmp_fm_struct)
     791        1312 :          CALL cp_fm_struct_release(tmp_fm_struct)
     792             :          CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
     793        1312 :                             mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
     794        1312 :          CALL choose_eigv_solver(overlap, evecs, evals)
     795       15160 :          mo_mag_min = MIN(MINVAL(evals), mo_mag_min)
     796       15160 :          mo_mag_max = MAX(MAXVAL(evals), mo_mag_max)
     797        1312 :          CALL cp_fm_release(overlap)
     798        1312 :          CALL cp_fm_release(evecs)
     799        5248 :          DEALLOCATE (evals)
     800             :       END DO
     801        1312 :       CALL timestop(handle)
     802             : 
     803        2624 :    END SUBROUTINE calculate_magnitude
     804             : 
     805             : ! **************************************************************************************************
     806             : !> \brief  Calculate KS eigenvalues starting from  OF MOS
     807             : !> \param mos ...
     808             : !> \param nspins ...
     809             : !> \param ks_rmpv ...
     810             : !> \param scf_control ...
     811             : !> \param mo_derivs ...
     812             : !> \param admm_env ...
     813             : !> \param hairy_probes ...
     814             : !> \param probe ...
     815             : !> \par History
     816             : !>         02.2013 moved from qs_scf_post_gpw
     817             : !>
     818             : ! **************************************************************************************************
     819         196 :    SUBROUTINE make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env, &
     820             :                           hairy_probes, probe)
     821             : 
     822             :       TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT)     :: mos
     823             :       INTEGER, INTENT(IN)                                :: nspins
     824             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_rmpv
     825             :       TYPE(scf_control_type), POINTER                    :: scf_control
     826             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mo_derivs
     827             :       TYPE(admm_type), OPTIONAL, POINTER                 :: admm_env
     828             :       LOGICAL, OPTIONAL                                  :: hairy_probes
     829             :       TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
     830             :          POINTER                                         :: probe
     831             : 
     832             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'make_mo_eig'
     833             : 
     834             :       INTEGER                                            :: handle, homo, ispin, nmo, output_unit
     835         196 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues
     836             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     837             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff_deriv
     838             : 
     839         196 :       CALL timeset(routineN, handle)
     840             : 
     841         196 :       NULLIFY (mo_coeff_deriv, mo_coeff, mo_eigenvalues)
     842         196 :       output_unit = cp_logger_get_default_io_unit()
     843             : 
     844         430 :       DO ispin = 1, nspins
     845             :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
     846         234 :                          eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
     847         234 :          IF (output_unit > 0) WRITE (output_unit, *) " "
     848         234 :          IF (output_unit > 0) WRITE (output_unit, *) " Eigenvalues of the occupied subspace spin ", ispin
     849             :          !      IF (homo .NE. nmo) THEN
     850             :          !         IF (output_unit>0) WRITE(output_unit,*)" and ",nmo-homo," added MO eigenvalues"
     851             :          !      END IF
     852         234 :          IF (output_unit > 0) WRITE (output_unit, *) "---------------------------------------------"
     853             : 
     854         234 :          IF (scf_control%use_ot) THEN
     855             :             ! Also rotate the OT derivs, since they are needed for force calculations
     856         128 :             IF (ASSOCIATED(mo_derivs)) THEN
     857         128 :                mo_coeff_deriv => mo_derivs(ispin)%matrix
     858             :             ELSE
     859           0 :                mo_coeff_deriv => NULL()
     860             :             END IF
     861             : 
     862             :             ! ** If we do ADMM, we add have to modify the kohn-sham matrix
     863         128 :             IF (PRESENT(admm_env)) THEN
     864           0 :                CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
     865             :             END IF
     866             : 
     867             :             CALL calculate_subspace_eigenvalues(mo_coeff, ks_rmpv(ispin)%matrix, mo_eigenvalues, &
     868             :                                                 scr=output_unit, ionode=output_unit > 0, do_rotation=.TRUE., &
     869         128 :                                                 co_rotate_dbcsr=mo_coeff_deriv)
     870             : 
     871             :             ! ** If we do ADMM, we restore the original kohn-sham matrix
     872         128 :             IF (PRESENT(admm_env)) THEN
     873           0 :                CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
     874             :             END IF
     875             :          ELSE
     876         106 :             IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(1:homo)
     877             :          END IF
     878         234 :          IF (.NOT. scf_control%diagonalization%mom) THEN
     879         234 :             IF (PRESENT(hairy_probes) .EQV. .TRUE.) THEN
     880           0 :                scf_control%smear%do_smear = .FALSE.
     881             :                CALL set_mo_occupation(mo_set=mos(ispin), &
     882             :                                       smear=scf_control%smear, &
     883           0 :                                       probe=probe)
     884             :             ELSE
     885         234 :                CALL set_mo_occupation(mo_set=mos(ispin), smear=scf_control%smear)
     886             :             END IF
     887             :          END IF
     888         234 :          IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') &
     889         547 :             "Fermi Energy [eV] :", mos(ispin)%mu*evolt
     890             :       END DO
     891             : 
     892         196 :       CALL timestop(handle)
     893             : 
     894         196 :    END SUBROUTINE make_mo_eig
     895             : 
     896             : END MODULE qs_mo_methods

Generated by: LCOV version 1.15