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