1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief lower level routines for linear scaling SCF 8!> \par History 9!> 2010.10 created [Joost VandeVondele] 10!> \author Joost VandeVondele 11! ************************************************************************************************** 12MODULE dm_ls_scf_methods 13 USE arnoldi_api, ONLY: arnoldi_extremal 14 USE cp_log_handling, ONLY: cp_get_default_logger,& 15 cp_logger_get_default_unit_nr,& 16 cp_logger_type 17 USE dbcsr_api, ONLY: & 18 dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_dot, & 19 dbcsr_filter, dbcsr_finalize, dbcsr_frobenius_norm, dbcsr_get_data_type, & 20 dbcsr_get_occupation, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, & 21 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, & 22 dbcsr_put_block, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_trace, dbcsr_type, & 23 dbcsr_type_no_symmetry, dbcsr_type_real_4 24 USE dm_ls_scf_qs, ONLY: matrix_qs_to_ls 25 USE dm_ls_scf_types, ONLY: ls_cluster_atomic,& 26 ls_mstruct_type,& 27 ls_scf_env_type 28 USE input_constants, ONLY: & 29 ls_cluster_atomic, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, & 30 ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_sign_ns, & 31 ls_scf_sign_proot, ls_scf_sign_submatrix, ls_scf_submatrix_sign_ns 32 USE iterate_matrix, ONLY: invert_Hotelling,& 33 matrix_sign_Newton_Schulz,& 34 matrix_sign_proot,& 35 matrix_sign_submatrix,& 36 matrix_sqrt_Newton_Schulz,& 37 matrix_sqrt_proot 38 USE kinds, ONLY: dp,& 39 int_8,& 40 sp 41 USE machine, ONLY: m_flush,& 42 m_walltime 43 USE mathlib, ONLY: abnormal_value 44#include "./base/base_uses.f90" 45 46 IMPLICIT NONE 47 48 PRIVATE 49 50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_methods' 51 52 PUBLIC :: ls_scf_init_matrix_S 53 PUBLIC :: density_matrix_sign, density_matrix_sign_fixed_mu 54 PUBLIC :: apply_matrix_preconditioner 55 PUBLIC :: density_matrix_trs4, density_matrix_tc2, compute_homo_lumo 56 57CONTAINS 58 59! ************************************************************************************************** 60!> \brief initialize S matrix related properties (sqrt, inverse...) 61!> Might be factored-out since this seems common code with the other SCF. 62!> \param matrix_s ... 63!> \param ls_scf_env ... 64!> \par History 65!> 2010.10 created [Joost VandeVondele] 66!> \author Joost VandeVondele 67! ************************************************************************************************** 68 SUBROUTINE ls_scf_init_matrix_S(matrix_s, ls_scf_env) 69 TYPE(dbcsr_type) :: matrix_s 70 TYPE(ls_scf_env_type) :: ls_scf_env 71 72 CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_init_matrix_S', & 73 routineP = moduleN//':'//routineN 74 75 INTEGER :: handle, unit_nr 76 REAL(KIND=dp) :: frob_matrix, frob_matrix_base 77 TYPE(cp_logger_type), POINTER :: logger 78 TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2 79 80 CALL timeset(routineN, handle) 81 82 ! get a useful output_unit 83 logger => cp_get_default_logger() 84 IF (logger%para_env%ionode) THEN 85 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 86 ELSE 87 unit_nr = -1 88 ENDIF 89 90 ! make our own copy of S 91 IF (ls_scf_env%has_unit_metric) THEN 92 CALL dbcsr_set(ls_scf_env%matrix_s, 0.0_dp) 93 CALL dbcsr_add_on_diag(ls_scf_env%matrix_s, 1.0_dp) 94 ELSE 95 CALL matrix_qs_to_ls(ls_scf_env%matrix_s, matrix_s, ls_scf_env%ls_mstruct, covariant=.TRUE.) 96 ENDIF 97 98 CALL dbcsr_filter(ls_scf_env%matrix_s, ls_scf_env%eps_filter) 99 100 ! needs a preconditioner for S 101 IF (ls_scf_env%has_s_preconditioner) THEN 102 CALL dbcsr_create(ls_scf_env%matrix_bs_sqrt, template=ls_scf_env%matrix_s, & 103 matrix_type=dbcsr_type_no_symmetry) 104 CALL dbcsr_create(ls_scf_env%matrix_bs_sqrt_inv, template=ls_scf_env%matrix_s, & 105 matrix_type=dbcsr_type_no_symmetry) 106 CALL compute_matrix_preconditioner(ls_scf_env%matrix_s, & 107 ls_scf_env%s_preconditioner_type, ls_scf_env%ls_mstruct, & 108 ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv, & 109 ls_scf_env%eps_filter, ls_scf_env%s_sqrt_order, & 110 ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos) 111 ENDIF 112 113 ! precondition S 114 IF (ls_scf_env%has_s_preconditioner) THEN 115 CALL apply_matrix_preconditioner(ls_scf_env%matrix_s, "forward", & 116 ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv) 117 ENDIF 118 119 ! compute sqrt(S) and inv(sqrt(S)) 120 IF (ls_scf_env%use_s_sqrt) THEN 121 122 CALL dbcsr_create(ls_scf_env%matrix_s_sqrt, template=ls_scf_env%matrix_s, & 123 matrix_type=dbcsr_type_no_symmetry) 124 CALL dbcsr_create(ls_scf_env%matrix_s_sqrt_inv, template=ls_scf_env%matrix_s, & 125 matrix_type=dbcsr_type_no_symmetry) 126 127 SELECT CASE (ls_scf_env%s_sqrt_method) 128 CASE (ls_s_sqrt_proot) 129 CALL matrix_sqrt_proot(ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_s_sqrt_inv, & 130 ls_scf_env%matrix_s, ls_scf_env%eps_filter, & 131 ls_scf_env%s_sqrt_order, & 132 ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos, & 133 symmetrize=.TRUE.) 134 CASE (ls_s_sqrt_ns) 135 CALL matrix_sqrt_Newton_Schulz(ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_s_sqrt_inv, & 136 ls_scf_env%matrix_s, ls_scf_env%eps_filter, & 137 ls_scf_env%s_sqrt_order, & 138 ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos) 139 CASE DEFAULT 140 CPABORT("Unknown sqrt method.") 141 END SELECT 142 143 IF (ls_scf_env%check_s_inv) THEN 144 CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, & 145 matrix_type=dbcsr_type_no_symmetry) 146 CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, & 147 matrix_type=dbcsr_type_no_symmetry) 148 149 CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_s, & 150 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter) 151 152 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, & 153 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter) 154 155 frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp2) 156 CALL dbcsr_add_on_diag(matrix_tmp2, -1.0_dp) 157 frob_matrix = dbcsr_frobenius_norm(matrix_tmp2) 158 IF (unit_nr > 0) THEN 159 WRITE (unit_nr, *) "Error for (inv(sqrt(S))*S*inv(sqrt(S))-I)", frob_matrix/frob_matrix_base 160 ENDIF 161 162 CALL dbcsr_release(matrix_tmp1) 163 CALL dbcsr_release(matrix_tmp2) 164 ENDIF 165 ENDIF 166 167 ! compute the inverse of S 168 IF (ls_scf_env%needs_s_inv) THEN 169 CALL dbcsr_create(ls_scf_env%matrix_s_inv, template=ls_scf_env%matrix_s, & 170 matrix_type=dbcsr_type_no_symmetry) 171 IF (.NOT. ls_scf_env%use_s_sqrt) THEN 172 CALL invert_Hotelling(ls_scf_env%matrix_s_inv, ls_scf_env%matrix_s, ls_scf_env%eps_filter) 173 ELSE 174 CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_s_sqrt_inv, & 175 0.0_dp, ls_scf_env%matrix_s_inv, filter_eps=ls_scf_env%eps_filter) 176 ENDIF 177 IF (ls_scf_env%check_s_inv) THEN 178 CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, & 179 matrix_type=dbcsr_type_no_symmetry) 180 CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_inv, ls_scf_env%matrix_s, & 181 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter) 182 frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp1) 183 CALL dbcsr_add_on_diag(matrix_tmp1, -1.0_dp) 184 frob_matrix = dbcsr_frobenius_norm(matrix_tmp1) 185 IF (unit_nr > 0) THEN 186 WRITE (unit_nr, *) "Error for (inv(S)*S-I)", frob_matrix/frob_matrix_base 187 ENDIF 188 CALL dbcsr_release(matrix_tmp1) 189 ENDIF 190 ENDIF 191 192 CALL timestop(handle) 193 END SUBROUTINE ls_scf_init_matrix_s 194 195! ************************************************************************************************** 196!> \brief compute for a block positive definite matrix s (bs) 197!> the sqrt(bs) and inv(sqrt(bs)) 198!> \param matrix_s ... 199!> \param preconditioner_type ... 200!> \param ls_mstruct ... 201!> \param matrix_bs_sqrt ... 202!> \param matrix_bs_sqrt_inv ... 203!> \param threshold ... 204!> \param order ... 205!> \param eps_lanczos ... 206!> \param max_iter_lanczos ... 207!> \par History 208!> 2010.10 created [Joost VandeVondele] 209!> \author Joost VandeVondele 210! ************************************************************************************************** 211 SUBROUTINE compute_matrix_preconditioner(matrix_s, preconditioner_type, ls_mstruct, & 212 matrix_bs_sqrt, matrix_bs_sqrt_inv, threshold, order, eps_lanczos, max_iter_lanczos) 213 214 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s 215 INTEGER :: preconditioner_type 216 TYPE(ls_mstruct_type) :: ls_mstruct 217 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_bs_sqrt, matrix_bs_sqrt_inv 218 REAL(KIND=dp) :: threshold 219 INTEGER :: order 220 REAL(KIND=dp) :: eps_lanczos 221 INTEGER :: max_iter_lanczos 222 223 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_matrix_preconditioner', & 224 routineP = moduleN//':'//routineN 225 226 INTEGER :: datatype, handle, iblock_col, iblock_row 227 LOGICAL :: block_needed 228 REAL(dp), DIMENSION(:, :), POINTER :: block_dp 229 REAL(sp), DIMENSION(:, :), POINTER :: block_sp 230 TYPE(dbcsr_iterator_type) :: iter 231 TYPE(dbcsr_type) :: matrix_bs 232 233 CALL timeset(routineN, handle) 234 235 datatype = dbcsr_get_data_type(matrix_s) ! could be single or double precision 236 237 ! first generate a block diagonal copy of s 238 CALL dbcsr_create(matrix_bs, template=matrix_s) 239 240 SELECT CASE (preconditioner_type) 241 CASE (ls_s_preconditioner_none) 242 CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular) 243 CALL dbcsr_iterator_start(iter, matrix_s) 244 DO WHILE (dbcsr_iterator_blocks_left(iter)) 245 IF (datatype == dbcsr_type_real_4) THEN 246 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, block_sp) 247 ELSE 248 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, block_dp) 249 ENDIF 250 251 ! do we need the block ? 252 ! this depends on the preconditioner, but also the matrix clustering method employed 253 ! for a clustered matrix, right now, we assume that atomic and molecular preconditioners 254 ! are actually the same, and only require that the diagonal blocks (clustered) are present 255 256 block_needed = .FALSE. 257 258 IF (iblock_row == iblock_col) THEN 259 block_needed = .TRUE. 260 ELSE 261 IF (preconditioner_type == ls_s_preconditioner_molecular .AND. & 262 ls_mstruct%cluster_type == ls_cluster_atomic) THEN 263 IF (ls_mstruct%atom_to_molecule(iblock_row) == ls_mstruct%atom_to_molecule(iblock_col)) block_needed = .TRUE. 264 ENDIF 265 ENDIF 266 267 ! add it 268 IF (block_needed) THEN 269 IF (datatype == dbcsr_type_real_4) THEN 270 CALL dbcsr_put_block(matrix=matrix_bs, row=iblock_row, col=iblock_col, block=block_sp) 271 ELSE 272 CALL dbcsr_put_block(matrix=matrix_bs, row=iblock_row, col=iblock_col, block=block_dp) 273 ENDIF 274 ENDIF 275 276 ENDDO 277 CALL dbcsr_iterator_stop(iter) 278 END SELECT 279 280 CALL dbcsr_finalize(matrix_bs) 281 282 SELECT CASE (preconditioner_type) 283 CASE (ls_s_preconditioner_none) 284 ! for now make it a simple identity matrix 285 CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs) 286 CALL dbcsr_set(matrix_bs_sqrt, 0.0_dp) 287 CALL dbcsr_add_on_diag(matrix_bs_sqrt, 1.0_dp) 288 289 ! for now make it a simple identity matrix 290 CALL dbcsr_copy(matrix_bs_sqrt_inv, matrix_bs) 291 CALL dbcsr_set(matrix_bs_sqrt_inv, 0.0_dp) 292 CALL dbcsr_add_on_diag(matrix_bs_sqrt_inv, 1.0_dp) 293 CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular) 294 CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs) 295 CALL dbcsr_copy(matrix_bs_sqrt_inv, matrix_bs) 296 ! XXXXXXXXXXX 297 ! XXXXXXXXXXX the threshold here could be done differently, 298 ! XXXXXXXXXXX using eps_filter is reducing accuracy for no good reason, this is cheap 299 ! XXXXXXXXXXX 300 CALL matrix_sqrt_Newton_Schulz(matrix_bs_sqrt, matrix_bs_sqrt_inv, matrix_bs, & 301 threshold=MIN(threshold, 1.0E-10_dp), order=order, & 302 eps_lanczos=eps_lanczos, max_iter_lanczos=max_iter_lanczos) 303 END SELECT 304 305 CALL dbcsr_release(matrix_bs) 306 307 CALL timestop(handle) 308 309 END SUBROUTINE compute_matrix_preconditioner 310 311! ************************************************************************************************** 312!> \brief apply a preconditioner either 313!> forward (precondition) inv(sqrt(bs)) * A * inv(sqrt(bs)) 314!> backward (restore to old form) sqrt(bs) * A * sqrt(bs) 315!> \param matrix ... 316!> \param direction ... 317!> \param matrix_bs_sqrt ... 318!> \param matrix_bs_sqrt_inv ... 319!> \par History 320!> 2010.10 created [Joost VandeVondele] 321!> \author Joost VandeVondele 322! ************************************************************************************************** 323 SUBROUTINE apply_matrix_preconditioner(matrix, direction, matrix_bs_sqrt, matrix_bs_sqrt_inv) 324 325 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 326 CHARACTER(LEN=*) :: direction 327 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_bs_sqrt, matrix_bs_sqrt_inv 328 329 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_matrix_preconditioner', & 330 routineP = moduleN//':'//routineN 331 332 INTEGER :: handle 333 TYPE(dbcsr_type) :: matrix_tmp 334 335 CALL timeset(routineN, handle) 336 CALL dbcsr_create(matrix_tmp, template=matrix, matrix_type=dbcsr_type_no_symmetry) 337 338 SELECT CASE (direction) 339 CASE ("forward") 340 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt_inv, & 341 0.0_dp, matrix_tmp) 342 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt_inv, matrix_tmp, & 343 0.0_dp, matrix) 344 CASE ("backward") 345 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt, & 346 0.0_dp, matrix_tmp) 347 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt, matrix_tmp, & 348 0.0_dp, matrix) 349 CASE DEFAULT 350 CPABORT("") 351 END SELECT 352 353 CALL dbcsr_release(matrix_tmp) 354 355 CALL timestop(handle) 356 357 END SUBROUTINE apply_matrix_preconditioner 358 359! ************************************************************************************************** 360!> \brief compute the density matrix with a trace that is close to nelectron. 361!> take a mu as input, and improve by bisection as needed. 362!> \param matrix_p ... 363!> \param mu ... 364!> \param fixed_mu ... 365!> \param sign_method ... 366!> \param sign_order ... 367!> \param matrix_ks ... 368!> \param matrix_s ... 369!> \param matrix_s_inv ... 370!> \param nelectron ... 371!> \param threshold ... 372!> \param sign_symmetric ... 373!> \param submatrix_sign_method ... 374!> \param matrix_s_sqrt_inv ... 375!> \par History 376!> 2010.10 created [Joost VandeVondele] 377!> \author Joost VandeVondele 378! ************************************************************************************************** 379 SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, & 380 matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric, submatrix_sign_method, & 381 matrix_s_sqrt_inv) 382 383 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p 384 REAL(KIND=dp), INTENT(INOUT) :: mu 385 LOGICAL :: fixed_mu 386 INTEGER :: sign_method, sign_order 387 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s, matrix_s_inv 388 INTEGER, INTENT(IN) :: nelectron 389 REAL(KIND=dp), INTENT(IN) :: threshold 390 LOGICAL, OPTIONAL :: sign_symmetric 391 INTEGER, OPTIONAL :: submatrix_sign_method 392 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_s_sqrt_inv 393 394 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_sign', & 395 routineP = moduleN//':'//routineN 396 REAL(KIND=dp), PARAMETER :: initial_increment = 0.01_dp 397 398 INTEGER :: handle, iter, unit_nr, & 399 used_submatrix_sign_method 400 LOGICAL :: do_sign_symmetric, has_mu_high, & 401 has_mu_low 402 REAL(KIND=dp) :: increment, mu_high, mu_low, trace 403 TYPE(cp_logger_type), POINTER :: logger 404 405 CALL timeset(routineN, handle) 406 407 logger => cp_get_default_logger() 408 IF (logger%para_env%ionode) THEN 409 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 410 ELSE 411 unit_nr = -1 412 ENDIF 413 414 do_sign_symmetric = .FALSE. 415 IF (PRESENT(sign_symmetric)) do_sign_symmetric = sign_symmetric 416 417 used_submatrix_sign_method = ls_scf_submatrix_sign_ns 418 IF (PRESENT(submatrix_sign_method)) used_submatrix_sign_method = submatrix_sign_method 419 420 increment = initial_increment 421 422 has_mu_low = .FALSE. 423 has_mu_high = .FALSE. 424 425 ! bisect if both bounds are known, otherwise find the bounds with a linear search 426 DO iter = 1, 30 427 IF (has_mu_low .AND. has_mu_high) THEN 428 mu = (mu_low + mu_high)/2 429 IF (ABS(mu_high - mu_low) < threshold) EXIT 430 ENDIF 431 432 CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, & 433 matrix_ks, matrix_s, matrix_s_inv, threshold, & 434 do_sign_symmetric, used_submatrix_sign_method, & 435 matrix_s_sqrt_inv) 436 IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,I2,1X,F13.9,1X,F15.9)') & 437 "Density matrix: iter, mu, trace error: ", iter, mu, trace - nelectron 438 439 ! OK, we can skip early if we are as close as possible to the exact result 440 ! smaller differences should be considered 'noise' 441 IF (ABS(trace - nelectron) < 0.5_dp .OR. fixed_mu) EXIT 442 443 IF (trace < nelectron) THEN 444 mu_low = mu 445 mu = mu + increment 446 has_mu_low = .TRUE. 447 increment = increment*2 448 ELSE 449 mu_high = mu 450 mu = mu - increment 451 has_mu_high = .TRUE. 452 increment = increment*2 453 ENDIF 454 ENDDO 455 456 CALL timestop(handle) 457 458 END SUBROUTINE density_matrix_sign 459 460! ************************************************************************************************** 461!> \brief for a fixed mu, compute the corresponding density matrix and its trace 462!> \param matrix_p ... 463!> \param trace ... 464!> \param mu ... 465!> \param sign_method ... 466!> \param sign_order ... 467!> \param matrix_ks ... 468!> \param matrix_s ... 469!> \param matrix_s_inv ... 470!> \param threshold ... 471!> \param sign_symmetric ... 472!> \param submatrix_sign_method ... 473!> \param matrix_s_sqrt_inv ... 474!> \par History 475!> 2010.10 created [Joost VandeVondele] 476!> \author Joost VandeVondele 477! ************************************************************************************************** 478 SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, matrix_ks, & 479 matrix_s, matrix_s_inv, threshold, sign_symmetric, submatrix_sign_method, & 480 matrix_s_sqrt_inv) 481 482 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p 483 REAL(KIND=dp), INTENT(OUT) :: trace 484 REAL(KIND=dp), INTENT(INOUT) :: mu 485 INTEGER :: sign_method, sign_order 486 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s, matrix_s_inv 487 REAL(KIND=dp), INTENT(IN) :: threshold 488 LOGICAL :: sign_symmetric 489 INTEGER :: submatrix_sign_method 490 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_s_sqrt_inv 491 492 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_sign_fixed_mu', & 493 routineP = moduleN//':'//routineN 494 495 INTEGER :: handle, unit_nr 496 REAL(KIND=dp) :: frob_matrix 497 TYPE(cp_logger_type), POINTER :: logger 498 TYPE(dbcsr_type) :: matrix_p_ud, matrix_sign, matrix_sinv_ks, matrix_ssqrtinv_ks_ssqrtinv, & 499 matrix_ssqrtinv_ks_ssqrtinv2, matrix_tmp 500 501 CALL timeset(routineN, handle) 502 503 logger => cp_get_default_logger() 504 IF (logger%para_env%ionode) THEN 505 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 506 ELSE 507 unit_nr = -1 508 ENDIF 509 510 CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry) 511 512 IF (sign_symmetric) THEN 513 514 IF (.NOT. PRESENT(matrix_s_sqrt_inv)) & 515 CPABORT("Argument matrix_s_sqrt_inv required if sign_symmetric is set") 516 517 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry) 518 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry) 519 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, & 520 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold) 521 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, & 522 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold) 523 CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu) 524 525 SELECT CASE (sign_method) 526 CASE (ls_scf_sign_ns) 527 CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order) 528 CASE (ls_scf_sign_proot) 529 CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order) 530 CASE (ls_scf_sign_submatrix) 531 CALL matrix_sign_submatrix(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, submatrix_sign_method) 532 CASE DEFAULT 533 CPABORT("Unkown sign method.") 534 END SELECT 535 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv) 536 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2) 537 538 ELSE ! .NOT. sign_symmetric 539 ! get inv(S)*H-I*mu 540 CALL dbcsr_create(matrix_sinv_ks, template=matrix_s, matrix_type=dbcsr_type_no_symmetry) 541 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_inv, matrix_ks, & 542 0.0_dp, matrix_sinv_ks, filter_eps=threshold) 543 CALL dbcsr_add_on_diag(matrix_sinv_ks, -mu) 544 545 ! compute sign(inv(S)*H-I*mu) 546 SELECT CASE (sign_method) 547 CASE (ls_scf_sign_ns) 548 CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order) 549 CASE (ls_scf_sign_proot) 550 CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order) 551 CASE (ls_scf_sign_submatrix) 552 CALL matrix_sign_submatrix(matrix_sign, matrix_sinv_ks, threshold, sign_order, submatrix_sign_method) 553 CASE DEFAULT 554 CPABORT("Unkown sign method.") 555 END SELECT 556 CALL dbcsr_release(matrix_sinv_ks) 557 ENDIF 558 559 ! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu)) 560 CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry) 561 CALL dbcsr_copy(matrix_p_ud, matrix_sign) 562 CALL dbcsr_scale(matrix_p_ud, -0.5_dp) 563 CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp) 564 CALL dbcsr_release(matrix_sign) 565 566 ! we now have PS, lets get its trace 567 CALL dbcsr_trace(matrix_p_ud, trace) 568 569 ! we can also check it is idempotent PS*PS=PS 570 CALL dbcsr_create(matrix_tmp, template=matrix_s, matrix_type=dbcsr_type_no_symmetry) 571 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_p_ud, & 572 0.0_dp, matrix_tmp, filter_eps=threshold) 573 CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp) 574 frob_matrix = dbcsr_frobenius_norm(matrix_tmp) 575 CALL dbcsr_release(matrix_tmp) 576 IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix 577 578 IF (sign_symmetric) THEN 579 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, & 580 0.0_dp, matrix_p, filter_eps=threshold) 581 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s_sqrt_inv, & 582 0.0_dp, matrix_p, filter_eps=threshold) 583 ELSE 584 585 ! get P=PS*inv(S) 586 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_s_inv, & 587 0.0_dp, matrix_p, filter_eps=threshold) 588 ENDIF 589 CALL dbcsr_release(matrix_p_ud) 590 591 CALL timestop(handle) 592 593 END SUBROUTINE density_matrix_sign_fixed_mu 594 595! ************************************************************************************************** 596!> \brief compute the density matrix using a trace-resetting algorithm 597!> \param matrix_p ... 598!> \param matrix_ks ... 599!> \param matrix_s_sqrt_inv ... 600!> \param nelectron ... 601!> \param threshold ... 602!> \param e_homo ... 603!> \param e_lumo ... 604!> \param e_mu ... 605!> \param dynamic_threshold ... 606!> \param matrix_ks_deviation ... 607!> \param max_iter_lanczos ... 608!> \param eps_lanczos ... 609!> \param converged ... 610!> \par History 611!> 2012.06 created [Florian Thoele] 612!> \author Florian Thoele 613! ************************************************************************************************** 614 SUBROUTINE density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, & 615 nelectron, threshold, e_homo, e_lumo, e_mu, & 616 dynamic_threshold, matrix_ks_deviation, & 617 max_iter_lanczos, eps_lanczos, converged) 618 619 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p 620 TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv 621 INTEGER, INTENT(IN) :: nelectron 622 REAL(KIND=dp), INTENT(IN) :: threshold 623 REAL(KIND=dp), INTENT(INOUT) :: e_homo, e_lumo, e_mu 624 LOGICAL, INTENT(IN), OPTIONAL :: dynamic_threshold 625 TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: matrix_ks_deviation 626 INTEGER, INTENT(IN) :: max_iter_lanczos 627 REAL(KIND=dp), INTENT(IN) :: eps_lanczos 628 LOGICAL, INTENT(OUT), OPTIONAL :: converged 629 630 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_trs4', & 631 routineP = moduleN//':'//routineN 632 INTEGER, PARAMETER :: max_iter = 100 633 REAL(KIND=dp), PARAMETER :: gamma_max = 6.0_dp, gamma_min = 0.0_dp 634 635 INTEGER :: branch, estimated_steps, handle, i, j, & 636 unit_nr 637 INTEGER(kind=int_8) :: flop1, flop2 638 LOGICAL :: arnoldi_converged, do_dyn_threshold 639 REAL(KIND=dp) :: current_threshold, delta_n, eps_max, eps_min, est_threshold, frob_id, & 640 frob_x, gam, homo, lumo, max_eig, max_threshold, maxdev, maxev, min_eig, minev, mmin, mu, & 641 mu_a, mu_b, mu_c, mu_fa, mu_fc, occ_matrix, scaled_homo_bound, scaled_lumo_bound, t1, t2, & 642 trace_fx, trace_gx, xi 643 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gamma_values 644 TYPE(cp_logger_type), POINTER :: logger 645 TYPE(dbcsr_type) :: matrix_k0, matrix_x, matrix_x_nosym, & 646 matrix_xidsq, matrix_xsq, tmp_gx 647 648 IF (nelectron == 0) THEN 649 CALL dbcsr_set(matrix_p, 0.0_dp) 650 RETURN 651 ENDIF 652 653 CALL timeset(routineN, handle) 654 655 logger => cp_get_default_logger() 656 IF (logger%para_env%ionode) THEN 657 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 658 ELSE 659 unit_nr = -1 660 ENDIF 661 662 do_dyn_threshold = .FALSE. 663 IF (PRESENT(dynamic_threshold)) do_dyn_threshold = dynamic_threshold 664 665 IF (PRESENT(converged)) converged = .FALSE. 666 667 ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2 668 CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type="S") 669 670 ! at some points the non-symmetric version of x is required 671 CALL dbcsr_create(matrix_x_nosym, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry) 672 673 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, & 674 0.0_dp, matrix_x_nosym, filter_eps=threshold) 675 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x_nosym, matrix_s_sqrt_inv, & 676 0.0_dp, matrix_x, filter_eps=threshold) 677 CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym) 678 679 CALL dbcsr_create(matrix_k0, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry) 680 CALL dbcsr_copy(matrix_k0, matrix_x_nosym) 681 682 ! compute the deviation in the mixed matrix, as seen in the ortho basis 683 IF (do_dyn_threshold) THEN 684 CPASSERT(PRESENT(matrix_ks_deviation)) 685 CALL dbcsr_add(matrix_ks_deviation, matrix_x_nosym, -1.0_dp, 1.0_dp) 686 CALL arnoldi_extremal(matrix_ks_deviation, maxev, minev, max_iter=max_iter_lanczos, threshold=eps_lanczos, & 687 converged=arnoldi_converged) 688 maxdev = MAX(ABS(maxev), ABS(minev)) 689 IF (unit_nr > 0) THEN 690 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", arnoldi_converged 691 WRITE (unit_nr, '(T6,A,1X,F12.5)') "change in mixed matrix: ", maxdev 692 WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo + maxdev 693 WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo - maxdev 694 WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo - maxdev) - (e_homo + maxdev)) > 0 695 ENDIF 696 ! save the old mixed matrix 697 CALL dbcsr_copy(matrix_ks_deviation, matrix_x_nosym) 698 699 ENDIF 700 701 ! get largest/smallest eigenvalues for scaling 702 CALL arnoldi_extremal(matrix_x_nosym, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, & 703 converged=arnoldi_converged) 704 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", & 705 min_eig, max_eig, " converged: ", arnoldi_converged 706 eps_max = max_eig 707 eps_min = min_eig 708 709 ! scale KS matrix 710 IF (eps_max == eps_min) THEN 711 CALL dbcsr_scale(matrix_x, 1.0_dp/eps_max) 712 ELSE 713 CALL dbcsr_add_on_diag(matrix_x, -eps_max) 714 CALL dbcsr_scale(matrix_x, -1.0_dp/(eps_max - eps_min)) 715 ENDIF 716 717 current_threshold = threshold 718 IF (do_dyn_threshold) THEN 719 ! scale bounds for HOMO/LUMO 720 scaled_homo_bound = (eps_max - (e_homo + maxdev))/(eps_max - eps_min) 721 scaled_lumo_bound = (eps_max - (e_lumo - maxdev))/(eps_max - eps_min) 722 ENDIF 723 724 CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type="S") 725 726 CALL dbcsr_create(matrix_xidsq, template=matrix_ks, matrix_type="S") 727 728 CALL dbcsr_create(tmp_gx, template=matrix_ks, matrix_type="S") 729 730 ALLOCATE (gamma_values(max_iter)) 731 732 DO i = 1, max_iter 733 t1 = m_walltime() 734 flop1 = 0; flop2 = 0 735 736 ! get X*X 737 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, & 738 0.0_dp, matrix_xsq, & 739 filter_eps=current_threshold, flop=flop1) 740 741 ! intermediate use matrix_xidsq to compute = X*X-X 742 CALL dbcsr_copy(matrix_xidsq, matrix_x) 743 CALL dbcsr_add(matrix_xidsq, matrix_xsq, -1.0_dp, 1.0_dp) 744 frob_id = dbcsr_frobenius_norm(matrix_xidsq) 745 frob_x = dbcsr_frobenius_norm(matrix_x) 746 747 ! xidsq = (1-X)*(1-X) 748 ! use (1-x)*(1-x) = 1 + x*x - 2*x 749 CALL dbcsr_copy(matrix_xidsq, matrix_x) 750 CALL dbcsr_add(matrix_xidsq, matrix_xsq, -2.0_dp, 1.0_dp) 751 CALL dbcsr_add_on_diag(matrix_xidsq, 1.0_dp) 752 753 ! tmp_gx = 4X-3X*X 754 CALL dbcsr_copy(tmp_gx, matrix_x) 755 CALL dbcsr_add(tmp_gx, matrix_xsq, 4.0_dp, -3.0_dp) 756 757 ! get gamma 758 ! Tr(F) = Tr(XX*tmp_gx) Tr(G) is equivalent 759 CALL dbcsr_dot(matrix_xsq, matrix_xidsq, trace_gx) 760 CALL dbcsr_dot(matrix_xsq, tmp_gx, trace_fx) 761 762 ! if converged, and gam becomes noisy, fix it to 3, which results in a final McWeeny step. 763 ! do this only if the electron count is reasonable. 764 ! maybe tune if the current criterion is not good enough 765 delta_n = nelectron - trace_fx 766 ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ... 767 IF (((frob_id*frob_id) < (threshold*frob_x*frob_x)) .AND. (ABS(delta_n) < 0.5_dp)) THEN 768 gam = 3.0_dp 769 ELSE IF (ABS(delta_n) < 1e-14_dp) THEN 770 gam = 0.0_dp ! rare case of perfect electron count 771 ELSE 772 ! make sure, we don't divide by zero, as soon as gam is outside the interval gam_min,gam_max, it doesn't matter 773 gam = delta_n/MAX(trace_gx, ABS(delta_n)/100) 774 ENDIF 775 gamma_values(i) = gam 776 777 IF (unit_nr > 0 .AND. .FALSE.) THEN 778 WRITE (unit_nr, *) "trace_fx", trace_fx, "trace_gx", trace_gx, "gam", gam, & 779 "frob_id", frob_id, "conv", ABS(frob_id/frob_x) 780 ENDIF 781 782 IF (do_dyn_threshold) THEN 783 ! quantities used for dynamic thresholding, when the estimated gap is larger than zero 784 xi = (scaled_homo_bound - scaled_lumo_bound) 785 IF (xi > 0.0_dp) THEN 786 mmin = 0.5*(scaled_homo_bound + scaled_lumo_bound) 787 max_threshold = ABS(1 - 2*mmin)*xi 788 789 scaled_homo_bound = evaluate_trs4_polynomial(scaled_homo_bound, gamma_values(i:), 1) 790 scaled_lumo_bound = evaluate_trs4_polynomial(scaled_lumo_bound, gamma_values(i:), 1) 791 estimated_steps = estimate_steps(scaled_homo_bound, scaled_lumo_bound, threshold) 792 793 est_threshold = (threshold/(estimated_steps + i + 1))*xi/(1 + threshold/(estimated_steps + i + 1)) 794 est_threshold = MIN(max_threshold, est_threshold) 795 IF (i > 1) est_threshold = MAX(est_threshold, 0.1_dp*current_threshold) 796 current_threshold = est_threshold 797 ELSE 798 current_threshold = threshold 799 ENDIF 800 ENDIF 801 802 IF (gam > gamma_max) THEN 803 ! Xn+1 = 2X-X*X 804 CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp) 805 CALL dbcsr_filter(matrix_x, current_threshold) 806 branch = 1 807 ELSE IF (gam < gamma_min) THEN 808 ! Xn+1 = X*X 809 CALL dbcsr_copy(matrix_x, matrix_xsq) 810 branch = 2 811 ELSE 812 ! Xn+1 = F(X) + gam*G(X) 813 CALL dbcsr_add(tmp_gx, matrix_xidsq, 1.0_dp, gam) 814 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, tmp_gx, & 815 0.0_dp, matrix_x, & 816 flop=flop2, filter_eps=current_threshold) 817 branch = 3 818 ENDIF 819 820 occ_matrix = dbcsr_get_occupation(matrix_x) 821 t2 = m_walltime() 822 IF (unit_nr > 0) THEN 823 WRITE (unit_nr, & 824 '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TRS4 it ", & 825 i, occ_matrix, ABS(trace_gx), t2 - t1, & 826 (flop1 + flop2)/(1.0E6_dp*MAX(t2 - t1, 0.001_dp)), current_threshold 827 CALL m_flush(unit_nr) 828 ENDIF 829 830 IF (abnormal_value(trace_gx)) & 831 CPABORT("trace_gx is an abnormal value (NaN/Inf).") 832 833 ! a branch of 1 or 2 appears to lead to a less accurate electron number count and premature exit 834 ! if it turns out this does not exit because we get stuck in branch 1/2 for a reason we need to refine further 835 ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ... 836 IF ((frob_id*frob_id) < (threshold*frob_x*frob_x) .AND. branch == 3 .AND. (ABS(delta_n) < 0.5_dp)) THEN 837 IF (PRESENT(converged)) converged = .TRUE. 838 EXIT 839 ENDIF 840 841 END DO 842 843 occ_matrix = dbcsr_get_occupation(matrix_x) 844 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,I3,1X,F10.8,E12.3)') 'Final TRS4 iteration ', i, occ_matrix, ABS(trace_gx) 845 846 ! free some memory 847 CALL dbcsr_release(tmp_gx) 848 CALL dbcsr_release(matrix_xsq) 849 CALL dbcsr_release(matrix_xidsq) 850 851 ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5 852 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, & 853 0.0_dp, matrix_x_nosym, filter_eps=threshold) 854 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_x_nosym, & 855 0.0_dp, matrix_p, filter_eps=threshold) 856 857 ! calculate the chemical potential by doing a bisection of fk(x0)-0.5, where fk is evaluated using the stored values for gamma 858 ! E. Rubensson et al., Chem Phys Lett 432, 2006, 591-594 859 mu_a = 0.0_dp; mu_b = 1.0_dp; 860 mu_fa = evaluate_trs4_polynomial(mu_a, gamma_values, i - 1) - 0.5_dp 861 DO j = 1, 40 862 mu_c = 0.5*(mu_a + mu_b) 863 mu_fc = evaluate_trs4_polynomial(mu_c, gamma_values, i - 1) - 0.5_dp ! i-1 because in the last iteration, only convergence is checked 864 IF (ABS(mu_fc) < 1.0E-6_dp .OR. (mu_b - mu_a)/2 < 1.0E-6_dp) EXIT !TODO: define threshold values 865 866 IF (mu_fc*mu_fa > 0) THEN 867 mu_a = mu_c 868 mu_fa = mu_fc 869 ELSE 870 mu_b = mu_c 871 ENDIF 872 ENDDO 873 mu = (eps_min - eps_max)*mu_c + eps_max 874 DEALLOCATE (gamma_values) 875 IF (unit_nr > 0) THEN 876 WRITE (unit_nr, '(T6,A,1X,F12.5)') 'Chemical potential (mu): ', mu 877 ENDIF 878 e_mu = mu 879 880 IF (do_dyn_threshold) THEN 881 CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym) 882 CALL compute_homo_lumo(matrix_k0, matrix_x_nosym, eps_min, eps_max, & 883 threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr) 884 e_homo = homo 885 e_lumo = lumo 886 ENDIF 887 888 CALL dbcsr_release(matrix_x) 889 CALL dbcsr_release(matrix_x_nosym) 890 CALL dbcsr_release(matrix_k0) 891 CALL timestop(handle) 892 893 END SUBROUTINE density_matrix_trs4 894 895! ************************************************************************************************** 896!> \brief compute the density matrix using a non monotonic trace conserving 897!> algorithm based on SIAM DOI. 10.1137/130911585. 898!> 2014.04 created [Jonathan Mullin] 899!> \param matrix_p ... 900!> \param matrix_ks ... 901!> \param matrix_s_sqrt_inv ... 902!> \param nelectron ... 903!> \param threshold ... 904!> \param e_homo ... 905!> \param e_lumo ... 906!> \param non_monotonic ... 907!> \param eps_lanczos ... 908!> \param max_iter_lanczos ... 909!> \author Jonathan Mullin 910! ************************************************************************************************** 911 SUBROUTINE density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, & 912 nelectron, threshold, e_homo, e_lumo, & 913 non_monotonic, eps_lanczos, max_iter_lanczos) 914 915 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p 916 TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv 917 INTEGER, INTENT(IN) :: nelectron 918 REAL(KIND=dp), INTENT(IN) :: threshold 919 REAL(KIND=dp), INTENT(INOUT) :: e_homo, e_lumo 920 LOGICAL, INTENT(IN), OPTIONAL :: non_monotonic 921 REAL(KIND=dp), INTENT(IN) :: eps_lanczos 922 INTEGER, INTENT(IN) :: max_iter_lanczos 923 924 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_tc2', & 925 routineP = moduleN//':'//routineN 926 INTEGER, PARAMETER :: max_iter = 100 927 928 INTEGER :: handle, i, j, k, unit_nr 929 INTEGER(kind=int_8) :: flop1, flop2 930 LOGICAL :: converged, do_non_monotonic 931 REAL(KIND=dp) :: beta, betaB, eps_max, eps_min, gama, & 932 max_eig, min_eig, occ_matrix, t1, t2, & 933 trace_fx, trace_gx 934 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, lambda, nu, poly, wu, X, Y 935 TYPE(cp_logger_type), POINTER :: logger 936 TYPE(dbcsr_type) :: matrix_tmp, matrix_x, matrix_xsq 937 938 CALL timeset(routineN, handle) 939 940 logger => cp_get_default_logger() 941 IF (logger%para_env%ionode) THEN 942 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 943 ELSE 944 unit_nr = -1 945 ENDIF 946 947 do_non_monotonic = .FALSE. 948 IF (PRESENT(non_monotonic)) do_non_monotonic = non_monotonic 949 950 ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2 951 CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry) 952 953 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, & 954 0.0_dp, matrix_x, filter_eps=threshold) 955 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, & 956 0.0_dp, matrix_x, filter_eps=threshold) 957 958 IF (unit_nr > 0) THEN 959 WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo 960 WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo 961 WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo) - (e_homo)) > 0 962 ENDIF 963 964 ! get largest/smallest eigenvalues for scaling 965 CALL arnoldi_extremal(matrix_x, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, & 966 converged=converged) 967 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", & 968 min_eig, max_eig, " converged: ", converged 969 970 eps_max = max_eig 971 eps_min = min_eig 972 973 ! scale KS matrix 974 CALL dbcsr_scale(matrix_x, -1.0_dp) 975 CALL dbcsr_add_on_diag(matrix_x, eps_max) 976 CALL dbcsr_scale(matrix_x, 1/(eps_max - eps_min)) 977 978 CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry) 979 CALL dbcsr_copy(matrix_xsq, matrix_x) 980 981 CALL dbcsr_create(matrix_tmp, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry) 982 983 ALLOCATE (poly(max_iter)) 984 ALLOCATE (nu(max_iter)) 985 ALLOCATE (wu(max_iter)) 986 ALLOCATE (alpha(max_iter)) 987 ALLOCATE (X(4)) 988 ALLOCATE (Y(4)) 989 ALLOCATE (lambda(4)) 990 991! Controls over the non_monotonic bounds, First if low gap, bias slightly 992 beta = (eps_max - ABS(e_lumo))/(eps_max - eps_min) 993 betaB = (eps_max + ABS(e_homo))/(eps_max - eps_min) 994 995 IF ((beta - betaB) < 0.005_dp) THEN 996 beta = beta - 0.002_dp 997 betaB = betaB + 0.002_dp 998 ENDIF 999! Check if input specifies to use monotonic bounds. 1000 IF (.NOT. do_non_monotonic) THEN 1001 beta = 0.0_dp 1002 betaB = 1.0_dp 1003 ENDIF 1004! initial SCF cycle has no reliable estimate of homo/lumo, force monotinic bounds. 1005 IF (e_homo == 0.0_dp) THEN 1006 beta = 0.0_dp 1007 BetaB = 1.0_dp 1008 ENDIF 1009 1010 ! init to take true branch first 1011 trace_fx = nelectron 1012 trace_gx = 0 1013 1014 DO i = 1, max_iter 1015 t1 = m_walltime() 1016 flop1 = 0; flop2 = 0 1017 1018 IF (ABS(trace_fx - nelectron) <= ABS(trace_gx - nelectron)) THEN 1019! Xn+1 = (aX+ (1-a)I)^2 1020 poly(i) = 1.0_dp 1021 alpha(i) = 2.0_dp/(2.0_dp - beta) 1022 1023 CALL dbcsr_scale(matrix_x, alpha(i)) 1024 CALL dbcsr_add_on_diag(matrix_x, 1.0_dp - alpha(i)) 1025 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, & 1026 0.0_dp, matrix_xsq, & 1027 filter_eps=threshold, flop=flop1) 1028 1029!save X for control variables 1030 CALL dbcsr_copy(matrix_tmp, matrix_x) 1031 1032 CALL dbcsr_copy(matrix_x, matrix_xsq) 1033 1034 beta = (1.0_dp - alpha(i)) + alpha(i)*beta 1035 beta = beta*beta 1036 betaB = (1.0_dp - alpha(i)) + alpha(i)*betaB 1037 betaB = betaB*betaB 1038 ELSE 1039! Xn+1 = 2aX-a^2*X^2 1040 poly(i) = 0.0_dp 1041 alpha(i) = 2.0_dp/(1.0_dp + betaB) 1042 1043 CALL dbcsr_scale(matrix_x, alpha(i)) 1044 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, & 1045 0.0_dp, matrix_xsq, & 1046 filter_eps=threshold, flop=flop1) 1047 1048!save X for control variables 1049 CALL dbcsr_copy(matrix_tmp, matrix_x) 1050! 1051 CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp) 1052 1053 beta = alpha(i)*beta 1054 beta = 2.0_dp*beta - beta*beta 1055 betaB = alpha(i)*betaB 1056 betaB = 2.0_dp*betaB - betaB*betaB 1057 1058 ENDIF 1059 occ_matrix = dbcsr_get_occupation(matrix_x) 1060 t2 = m_walltime() 1061 IF (unit_nr > 0) THEN 1062 WRITE (unit_nr, & 1063 '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TC2 it ", & 1064 i, occ_matrix, t2 - t1, & 1065 (flop1 + flop2)/(1.0E6_dp*(t2 - t1)), threshold 1066 CALL m_flush(unit_nr) 1067 ENDIF 1068 1069! calculate control terms 1070 CALL dbcsr_trace(matrix_xsq, trace_fx) 1071 1072! intermediate use matrix_xsq compute X- X*X , temorarily use trace_gx 1073 CALL dbcsr_add(matrix_xsq, matrix_tmp, -1.0_dp, 1.0_dp) 1074 CALL dbcsr_trace(matrix_xsq, trace_gx) 1075 nu(i) = dbcsr_frobenius_norm(matrix_xsq) 1076 wu(i) = trace_gx 1077 1078! intermediate use matrix_xsq to compute = 2X - X*X 1079 CALL dbcsr_add(matrix_xsq, matrix_tmp, 1.0_dp, 1.0_dp) 1080 CALL dbcsr_trace(matrix_xsq, trace_gx) 1081! TC2 has quadratic convergence, using the frobeniums norm as an idempotency deviation test. 1082 IF (ABS(nu(i)) < (threshold)) EXIT 1083 END DO 1084 1085 occ_matrix = dbcsr_get_occupation(matrix_x) 1086 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,I3,1X,1F10.8,1X,1F10.8)') 'Final TC2 iteration ', i, occ_matrix, ABS(nu(i)) 1087 1088 CALL dbcsr_release(matrix_xsq) 1089 CALL dbcsr_release(matrix_tmp) 1090 1091 ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5 1092 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, & 1093 0.0_dp, matrix_p, filter_eps=threshold) 1094 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p, & 1095 0.0_dp, matrix_p, filter_eps=threshold) 1096 1097 ! ALGO 3 from. SIAM DOI. 10.1137/130911585 1098 X(1) = 1.0_dp 1099 X(2) = 1.0_dp 1100 X(3) = 0.0_dp 1101 X(4) = 0.0_dp 1102 gama = 6.0_dp - 4.0_dp*(SQRT(2.0_dp)) 1103 gama = gama - gama*gama 1104 DO WHILE (nu(i) < gama) 1105 ! safeguard against negative root, is skipping correct? 1106 IF (wu(i) < 1.0e-14_dp) THEN 1107 i = i - 1 1108 CYCLE 1109 END IF 1110 IF ((1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)) < 0.0_dp) THEN 1111 i = i - 1 1112 CYCLE 1113 END IF 1114 Y(1) = 0.5_dp*(1.0_dp - SQRT(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i))) 1115 Y(2) = 0.5_dp*(1.0_dp - SQRT(1.0_dp - 4.0_dp*nu(i))) 1116 Y(3) = 0.5_dp*(1.0_dp + SQRT(1.0_dp - 4.0_dp*nu(i))) 1117 Y(4) = 0.5_dp*(1.0_dp + SQRT(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i))) 1118 Y(:) = MIN(1.0_dp, MAX(0.0_dp, Y(:))) 1119 DO j = i, 1, -1 1120 IF (poly(j) == 1.0_dp) THEN 1121 DO k = 1, 4 1122 Y(k) = SQRT(Y(k)) 1123 Y(k) = (Y(k) - 1.0_dp + alpha(j))/alpha(j) 1124 ENDDO ! end K 1125 ELSE 1126 DO k = 1, 4 1127 Y(k) = 1.0_dp - SQRT(1.0_dp - Y(k)) 1128 Y(k) = Y(k)/alpha(j) 1129 ENDDO ! end K 1130 ENDIF ! end poly 1131 ENDDO ! end j 1132 X(1) = MIN(X(1), Y(1)) 1133 X(2) = MIN(X(2), Y(2)) 1134 X(3) = MAX(X(3), Y(3)) 1135 X(4) = MAX(X(4), Y(4)) 1136 i = i - 1 1137 ENDDO ! end i 1138! lambda 1,2,3,4 are:: out lumo, in lumo, in homo, out homo 1139 DO k = 1, 4 1140 lambda(k) = eps_max - (eps_max - eps_min)*X(k) 1141 ENDDO ! end k 1142! END ALGO 3 from. SIAM DOI. 10.1137/130911585 1143 e_homo = lambda(4) 1144 e_lumo = lambda(1) 1145 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "outer homo/lumo/gap", e_homo, e_lumo, (e_lumo - e_homo) 1146 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "inner homo/lumo/gap", lambda(3), lambda(2), (lambda(2) - lambda(3)) 1147 1148 DEALLOCATE (poly) 1149 DEALLOCATE (nu) 1150 DEALLOCATE (wu) 1151 DEALLOCATE (alpha) 1152 DEALLOCATE (X) 1153 DEALLOCATE (Y) 1154 DEALLOCATE (lambda) 1155 1156 CALL dbcsr_release(matrix_x) 1157 CALL timestop(handle) 1158 1159 END SUBROUTINE density_matrix_tc2 1160 1161! ************************************************************************************************** 1162!> \brief compute the homo and lumo given a KS matrix and a density matrix in the orthonormalized basis 1163!> and the eps_min and eps_max, min and max eigenvalue of the ks matrix 1164!> \param matrix_k ... 1165!> \param matrix_p ... 1166!> \param eps_min ... 1167!> \param eps_max ... 1168!> \param threshold ... 1169!> \param max_iter_lanczos ... 1170!> \param eps_lanczos ... 1171!> \param homo ... 1172!> \param lumo ... 1173!> \param unit_nr ... 1174!> \par History 1175!> 2012.06 created [Florian Thoele] 1176!> \author Florian Thoele 1177! ************************************************************************************************** 1178 SUBROUTINE compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr) 1179 TYPE(dbcsr_type) :: matrix_k, matrix_p 1180 REAL(KIND=dp) :: eps_min, eps_max, threshold 1181 INTEGER, INTENT(IN) :: max_iter_lanczos 1182 REAL(KIND=dp), INTENT(IN) :: eps_lanczos 1183 REAL(KIND=dp) :: homo, lumo 1184 INTEGER :: unit_nr 1185 1186 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_homo_lumo', & 1187 routineP = moduleN//':'//routineN 1188 1189 LOGICAL :: converged 1190 REAL(KIND=dp) :: max_eig, min_eig, shift1, shift2 1191 TYPE(dbcsr_type) :: tmp1, tmp2, tmp3 1192 1193! temporary matrices used for HOMO/LUMO calculation 1194 1195 CALL dbcsr_create(tmp1, template=matrix_k, matrix_type=dbcsr_type_no_symmetry) 1196 1197 CALL dbcsr_create(tmp2, template=matrix_k, matrix_type=dbcsr_type_no_symmetry) 1198 1199 CALL dbcsr_create(tmp3, template=matrix_k, matrix_type=dbcsr_type_no_symmetry) 1200 1201 shift1 = -eps_min 1202 shift2 = eps_max 1203 1204 ! find largest ev of P*(K+shift*1), where shift is the neg. val. of the smallest ev of K 1205 CALL dbcsr_copy(tmp2, matrix_k) 1206 CALL dbcsr_add_on_diag(tmp2, shift1) 1207 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, tmp2, & 1208 0.0_dp, tmp1, filter_eps=threshold) 1209 CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, & 1210 threshold=eps_lanczos, max_iter=max_iter_lanczos) 1211 homo = max_eig - shift1 1212 IF (unit_nr > 0) THEN 1213 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged 1214 ENDIF 1215 1216 ! -(1-P)*(K-shift*1) = (1-P)*(shift*1 - K), where shift is the largest ev of K 1217 CALL dbcsr_copy(tmp3, matrix_p) 1218 CALL dbcsr_scale(tmp3, -1.0_dp) 1219 CALL dbcsr_add_on_diag(tmp3, 1.0_dp) !tmp3 = 1-P 1220 CALL dbcsr_copy(tmp2, matrix_k) 1221 CALL dbcsr_add_on_diag(tmp2, -shift2) 1222 CALL dbcsr_multiply("N", "N", -1.0_dp, tmp3, tmp2, & 1223 0.0_dp, tmp1, filter_eps=threshold) 1224 CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, & 1225 threshold=eps_lanczos, max_iter=max_iter_lanczos) 1226 lumo = -max_eig + shift2 1227 1228 IF (unit_nr > 0) THEN 1229 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged 1230 WRITE (unit_nr, '(T6,A,1X,3F12.5)') 'HOMO/LUMO/gap', homo, lumo, lumo - homo 1231 ENDIF 1232 CALL dbcsr_release(tmp1) 1233 CALL dbcsr_release(tmp2) 1234 CALL dbcsr_release(tmp3) 1235 1236 END SUBROUTINE compute_homo_lumo 1237 1238! ************************************************************************************************** 1239!> \brief ... 1240!> \param x ... 1241!> \param gamma_values ... 1242!> \param i ... 1243!> \return ... 1244! ************************************************************************************************** 1245 FUNCTION evaluate_trs4_polynomial(x, gamma_values, i) RESULT(xr) 1246 REAL(KIND=dp) :: x 1247 REAL(KIND=dp), DIMENSION(:) :: gamma_values 1248 INTEGER :: i 1249 REAL(KIND=dp) :: xr 1250 1251 REAL(KIND=dp), PARAMETER :: gam_max = 6.0_dp, gam_min = 0.0_dp 1252 1253 INTEGER :: k 1254 1255 xr = x 1256 DO k = 1, i 1257 IF (gamma_values(k) > gam_max) THEN 1258 xr = 2*xr - xr**2 1259 ELSE IF (gamma_values(k) < gam_min) THEN 1260 xr = xr**2 1261 ELSE 1262 xr = (xr*xr)*(4*xr - 3*xr*xr) + gamma_values(k)*xr*xr*((1 - xr)**2) 1263 ENDIF 1264 ENDDO 1265 END FUNCTION evaluate_trs4_polynomial 1266 1267! ************************************************************************************************** 1268!> \brief ... 1269!> \param homo ... 1270!> \param lumo ... 1271!> \param threshold ... 1272!> \return ... 1273! ************************************************************************************************** 1274 FUNCTION estimate_steps(homo, lumo, threshold) RESULT(steps) 1275 REAL(KIND=dp) :: homo, lumo, threshold 1276 INTEGER :: steps 1277 1278 INTEGER :: i 1279 REAL(KIND=dp) :: h, l, m 1280 1281 l = lumo 1282 h = homo 1283 1284 DO i = 1, 200 1285 IF (ABS(l) < threshold .AND. ABS(1 - h) < threshold) EXIT 1286 m = 0.5_dp*(h + l) 1287 IF (m > 0.5_dp) THEN 1288 h = h**2 1289 l = l**2 1290 ELSE 1291 h = 2*h - h**2 1292 l = 2*l - l**2 1293 ENDIF 1294 ENDDO 1295 steps = i - 1 1296 END FUNCTION estimate_steps 1297 1298END MODULE dm_ls_scf_methods 1299