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