1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Contains ADMM methods which only require the density matrix 8!> \par History 9!> 11.2014 created [Ole Schuett] 10!> \author Ole Schuett 11! ************************************************************************************************** 12MODULE admm_dm_methods 13 USE admm_dm_types, ONLY: admm_dm_type,& 14 mcweeny_history_type 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set 17 USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr 18 USE dbcsr_api, ONLY: & 19 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_frobenius_norm, dbcsr_get_block_p, & 20 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, & 21 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, & 22 dbcsr_scale, dbcsr_set, dbcsr_type 23 USE input_constants, ONLY: do_admm_basis_projection,& 24 do_admm_blocked_projection 25 USE iterate_matrix, ONLY: invert_Hotelling 26 USE kinds, ONLY: dp 27 USE pw_types, ONLY: pw_p_type 28 USE qs_collocate_density, ONLY: calculate_rho_elec 29 USE qs_ks_types, ONLY: get_ks_env,& 30 qs_ks_env_type 31 USE qs_rho_types, ONLY: qs_rho_get,& 32 qs_rho_set,& 33 qs_rho_type 34#include "./base/base_uses.f90" 35 36 IMPLICIT NONE 37 PRIVATE 38 39 PUBLIC :: admm_dm_calc_rho_aux, admm_dm_merge_ks_matrix 40 41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_dm_methods' 42 43CONTAINS 44 45! ************************************************************************************************** 46!> \brief Entry methods: Calculates auxilary density matrix from primary one. 47!> \param ks_env ... 48!> \author Ole Schuett 49! ************************************************************************************************** 50 SUBROUTINE admm_dm_calc_rho_aux(ks_env) 51 TYPE(qs_ks_env_type), POINTER :: ks_env 52 53 CHARACTER(len=*), PARAMETER :: routineN = 'admm_dm_calc_rho_aux', & 54 routineP = moduleN//':'//routineN 55 56 INTEGER :: handle 57 TYPE(admm_dm_type), POINTER :: admm_dm 58 59 NULLIFY (admm_dm) 60 CALL timeset(routineN, handle) 61 CALL get_ks_env(ks_env, admm_dm=admm_dm) 62 63 SELECT CASE (admm_dm%method) 64 CASE (do_admm_basis_projection) 65 CALL map_dm_projection(ks_env) 66 67 CASE (do_admm_blocked_projection) 68 CALL map_dm_blocked(ks_env) 69 70 CASE DEFAULT 71 CPABORT("admm_dm_calc_rho_aux: unknown method") 72 END SELECT 73 74 IF (admm_dm%purify) & 75 CALL purify_mcweeny(ks_env) 76 77 CALL update_rho_aux(ks_env) 78 79 CALL timestop(handle) 80 END SUBROUTINE admm_dm_calc_rho_aux 81 82! ************************************************************************************************** 83!> \brief Entry methods: Merges auxilary Kohn-Sham matrix into primary one. 84!> \param ks_env ... 85!> \author Ole Schuett 86! ************************************************************************************************** 87 SUBROUTINE admm_dm_merge_ks_matrix(ks_env) 88 TYPE(qs_ks_env_type), POINTER :: ks_env 89 90 CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_dm_merge_ks_matrix', & 91 routineP = moduleN//':'//routineN 92 93 INTEGER :: handle 94 TYPE(admm_dm_type), POINTER :: admm_dm 95 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge 96 97 CALL timeset(routineN, handle) 98 NULLIFY (admm_dm, matrix_ks_merge) 99 100 CALL get_ks_env(ks_env, admm_dm=admm_dm) 101 102 IF (admm_dm%purify) THEN 103 CALL revert_purify_mcweeny(ks_env, matrix_ks_merge) 104 ELSE 105 CALL get_ks_env(ks_env, matrix_ks_aux_fit=matrix_ks_merge) 106 ENDIF 107 108 SELECT CASE (admm_dm%method) 109 CASE (do_admm_basis_projection) 110 CALL merge_dm_projection(ks_env, matrix_ks_merge) 111 112 CASE (do_admm_blocked_projection) 113 CALL merge_dm_blocked(ks_env, matrix_ks_merge) 114 115 CASE DEFAULT 116 CPABORT("admm_dm_merge_ks_matrix: unknown method") 117 END SELECT 118 119 IF (admm_dm%purify) & 120 CALL dbcsr_deallocate_matrix_set(matrix_ks_merge) 121 122 CALL timestop(handle) 123 124 END SUBROUTINE admm_dm_merge_ks_matrix 125 126! ************************************************************************************************** 127!> \brief Calculates auxilary density matrix via basis projection. 128!> \param ks_env ... 129!> \author Ole Schuett 130! ************************************************************************************************** 131 SUBROUTINE map_dm_projection(ks_env) 132 TYPE(qs_ks_env_type), POINTER :: ks_env 133 134 CHARACTER(len=*), PARAMETER :: routineN = 'map_dm_projection', & 135 routineP = moduleN//':'//routineN 136 137 INTEGER :: ispin 138 LOGICAL :: s_mstruct_changed 139 REAL(KIND=dp) :: threshold 140 TYPE(admm_dm_type), POINTER :: admm_dm 141 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux, matrix_s_mixed, rho_ao, & 142 rho_ao_aux 143 TYPE(dbcsr_type) :: matrix_s_aux_inv, matrix_tmp 144 TYPE(dft_control_type), POINTER :: dft_control 145 TYPE(qs_rho_type), POINTER :: rho, rho_aux 146 147 NULLIFY (dft_control, admm_dm, matrix_s_aux, matrix_s_mixed, rho, rho_aux) 148 NULLIFY (rho_ao, rho_ao_aux) 149 150 CALL get_ks_env(ks_env, & 151 admm_dm=admm_dm, & 152 dft_control=dft_control, & 153 matrix_s_aux_fit=matrix_s_aux, & 154 matrix_s_aux_fit_vs_orb=matrix_s_mixed, & 155 s_mstruct_changed=s_mstruct_changed, & 156 rho=rho, & 157 rho_aux_fit=rho_aux) 158 159 CALL qs_rho_get(rho, rho_ao=rho_ao) 160 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux) 161 162 IF (s_mstruct_changed) THEN 163 ! Calculate A = S_aux^(-1) * S_mixed 164 CALL dbcsr_create(matrix_s_aux_inv, template=matrix_s_aux(1)%matrix, matrix_type="N") 165 threshold = MAX(admm_dm%eps_filter, 1.0e-12_dp) 166 CALL invert_Hotelling(matrix_s_aux_inv, matrix_s_aux(1)%matrix, threshold) 167 168 IF (.NOT. ASSOCIATED(admm_dm%matrix_A)) THEN 169 ALLOCATE (admm_dm%matrix_A) 170 CALL dbcsr_create(admm_dm%matrix_A, template=matrix_s_mixed(1)%matrix, matrix_type="N") 171 ENDIF 172 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_aux_inv, matrix_s_mixed(1)%matrix, & 173 0.0_dp, admm_dm%matrix_A) 174 CALL dbcsr_release(matrix_s_aux_inv) 175 ENDIF 176 177 ! Calculate P_aux = A * P * A^T 178 CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A) 179 DO ispin = 1, dft_control%nspins 180 CALL dbcsr_multiply("N", "N", 1.0_dp, admm_dm%matrix_A, rho_ao(ispin)%matrix, & 181 0.0_dp, matrix_tmp) 182 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, admm_dm%matrix_A, & 183 0.0_dp, rho_ao_aux(ispin)%matrix) 184 END DO 185 CALL dbcsr_release(matrix_tmp) 186 187 END SUBROUTINE map_dm_projection 188 189! ************************************************************************************************** 190!> \brief Calculates auxilary density matrix via blocking. 191!> \param ks_env ... 192!> \author Ole Schuett 193! ************************************************************************************************** 194 SUBROUTINE map_dm_blocked(ks_env) 195 TYPE(qs_ks_env_type), POINTER :: ks_env 196 197 INTEGER :: blk, iatom, ispin, jatom 198 LOGICAL :: found 199 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux 200 TYPE(admm_dm_type), POINTER :: admm_dm 201 TYPE(dbcsr_iterator_type) :: iter 202 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_aux 203 TYPE(dft_control_type), POINTER :: dft_control 204 TYPE(qs_rho_type), POINTER :: rho, rho_aux 205 206 NULLIFY (dft_control, admm_dm, rho, rho_aux, rho_ao, rho_ao_aux) 207 208 CALL get_ks_env(ks_env, & 209 admm_dm=admm_dm, & 210 dft_control=dft_control, & 211 rho=rho, & 212 rho_aux_fit=rho_aux) 213 214 CALL qs_rho_get(rho, rho_ao=rho_ao) 215 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux) 216 217 ! ** set blocked density matrix to 0 218 DO ispin = 1, dft_control%nspins 219 CALL dbcsr_set(rho_ao_aux(ispin)%matrix, 0.0_dp) 220 ! ** now loop through the list and copy corresponding blocks 221 CALL dbcsr_iterator_start(iter, rho_ao(ispin)%matrix) 222 DO WHILE (dbcsr_iterator_blocks_left(iter)) 223 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) 224 IF (admm_dm%block_map(iatom, jatom) == 1) THEN 225 CALL dbcsr_get_block_p(rho_ao_aux(ispin)%matrix, & 226 row=iatom, col=jatom, BLOCK=sparse_block_aux, found=found) 227 IF (found) & 228 sparse_block_aux = sparse_block 229 END IF 230 END DO 231 CALL dbcsr_iterator_stop(iter) 232 ENDDO 233 234 END SUBROUTINE map_dm_blocked 235 236! ************************************************************************************************** 237!> \brief Call calculate_rho_elec() for auxilary density 238!> \param ks_env ... 239! ************************************************************************************************** 240 SUBROUTINE update_rho_aux(ks_env) 241 TYPE(qs_ks_env_type), POINTER :: ks_env 242 243 INTEGER :: ispin 244 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_aux 245 TYPE(admm_dm_type), POINTER :: admm_dm 246 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux 247 TYPE(dft_control_type), POINTER :: dft_control 248 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g_aux, rho_r_aux 249 TYPE(qs_rho_type), POINTER :: rho_aux 250 251 NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux) 252 253 CALL get_ks_env(ks_env, & 254 admm_dm=admm_dm, & 255 dft_control=dft_control, & 256 rho_aux_fit=rho_aux) 257 258 CALL qs_rho_get(rho_aux, & 259 rho_ao=rho_ao_aux, & 260 rho_r=rho_r_aux, & 261 rho_g=rho_g_aux, & 262 tot_rho_r=tot_rho_r_aux) 263 264 DO ispin = 1, dft_control%nspins 265 CALL calculate_rho_elec(ks_env=ks_env, & 266 matrix_p=rho_ao_aux(ispin)%matrix, & 267 rho=rho_r_aux(ispin), & 268 rho_gspace=rho_g_aux(ispin), & 269 total_rho=tot_rho_r_aux(ispin), & 270 soft_valid=.FALSE., & 271 basis_type="AUX_FIT") 272 END DO 273 274 CALL qs_rho_set(rho_aux, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) 275 276 END SUBROUTINE update_rho_aux 277 278! ************************************************************************************************** 279!> \brief Merges auxilary Kohn-Sham matrix via basis projection. 280!> \param ks_env ... 281!> \param matrix_ks_merge Input: The KS matrix to be merged 282!> \author Ole Schuett 283! ************************************************************************************************** 284 SUBROUTINE merge_dm_projection(ks_env, matrix_ks_merge) 285 TYPE(qs_ks_env_type), POINTER :: ks_env 286 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge 287 288 INTEGER :: ispin 289 TYPE(admm_dm_type), POINTER :: admm_dm 290 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks 291 TYPE(dbcsr_type) :: matrix_tmp 292 TYPE(dft_control_type), POINTER :: dft_control 293 294 NULLIFY (admm_dm, dft_control, matrix_ks) 295 296 CALL get_ks_env(ks_env, & 297 admm_dm=admm_dm, & 298 dft_control=dft_control, & 299 matrix_ks=matrix_ks) 300 301 ! Calculate K += A^T * K_aux * A 302 CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A, matrix_type="N") 303 304 DO ispin = 1, dft_control%nspins 305 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ks_merge(ispin)%matrix, admm_dm%matrix_A, & 306 0.0_dp, matrix_tmp) 307 CALL dbcsr_multiply("T", "N", 1.0_dp, admm_dm%matrix_A, matrix_tmp, & 308 1.0_dp, matrix_ks(ispin)%matrix) 309 END DO 310 311 CALL dbcsr_release(matrix_tmp) 312 313 END SUBROUTINE merge_dm_projection 314 315! ************************************************************************************************** 316!> \brief Merges auxilary Kohn-Sham matrix via blocking. 317!> \param ks_env ... 318!> \param matrix_ks_merge Input: The KS matrix to be merged 319!> \author Ole Schuett 320! ************************************************************************************************** 321 SUBROUTINE merge_dm_blocked(ks_env, matrix_ks_merge) 322 TYPE(qs_ks_env_type), POINTER :: ks_env 323 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge 324 325 INTEGER :: blk, iatom, ispin, jatom 326 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block 327 TYPE(admm_dm_type), POINTER :: admm_dm 328 TYPE(dbcsr_iterator_type) :: iter 329 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks 330 TYPE(dft_control_type), POINTER :: dft_control 331 332 NULLIFY (admm_dm, dft_control, matrix_ks) 333 334 CALL get_ks_env(ks_env, & 335 admm_dm=admm_dm, & 336 dft_control=dft_control, & 337 matrix_ks=matrix_ks) 338 339 DO ispin = 1, dft_control%nspins 340 CALL dbcsr_iterator_start(iter, matrix_ks_merge(ispin)%matrix) 341 DO WHILE (dbcsr_iterator_blocks_left(iter)) 342 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) 343 IF (admm_dm%block_map(iatom, jatom) == 0) & 344 sparse_block = 0.0_dp 345 END DO 346 CALL dbcsr_iterator_stop(iter) 347 CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_merge(ispin)%matrix, 1.0_dp, 1.0_dp) 348 ENDDO 349 350 END SUBROUTINE merge_dm_blocked 351 352! ************************************************************************************************** 353!> \brief Apply McWeeny purification to auxilary density matrix 354!> \param ks_env ... 355!> \author Ole Schuett 356! ************************************************************************************************** 357 SUBROUTINE purify_mcweeny(ks_env) 358 TYPE(qs_ks_env_type), POINTER :: ks_env 359 360 CHARACTER(LEN=*), PARAMETER :: routineN = 'purify_mcweeny', routineP = moduleN//':'//routineN 361 362 INTEGER :: handle, ispin, istep, nspins, unit_nr 363 REAL(KIND=dp) :: frob_norm 364 TYPE(admm_dm_type), POINTER :: admm_dm 365 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, rho_ao_aux 366 TYPE(dbcsr_type) :: matrix_ps, matrix_psp, matrix_test 367 TYPE(dbcsr_type), POINTER :: matrix_p, matrix_s 368 TYPE(dft_control_type), POINTER :: dft_control 369 TYPE(mcweeny_history_type), POINTER :: history, new_hist_entry 370 TYPE(qs_rho_type), POINTER :: rho_aux_fit 371 372 CALL timeset(routineN, handle) 373 NULLIFY (dft_control, admm_dm, matrix_s_aux_fit, rho_aux_fit, new_hist_entry, & 374 matrix_p, matrix_s, rho_ao_aux) 375 376 unit_nr = cp_logger_get_default_unit_nr() 377 CALL get_ks_env(ks_env, & 378 dft_control=dft_control, & 379 admm_dm=admm_dm, & 380 matrix_s_aux_fit=matrix_s_aux_fit, & 381 rho_aux_fit=rho_aux_fit) 382 383 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux) 384 385 matrix_p => rho_ao_aux(1)%matrix 386 CALL dbcsr_create(matrix_PS, template=matrix_p, matrix_type="N") 387 CALL dbcsr_create(matrix_PSP, template=matrix_p, matrix_type="S") 388 CALL dbcsr_create(matrix_test, template=matrix_p, matrix_type="S") 389 390 nspins = dft_control%nspins 391 DO ispin = 1, nspins 392 matrix_p => rho_ao_aux(ispin)%matrix 393 matrix_s => matrix_s_aux_fit(1)%matrix 394 history => admm_dm%mcweeny_history(ispin)%p 395 IF (ASSOCIATED(history)) CPABORT("purify_dm_mcweeny: history already associated") 396 IF (nspins == 1) CALL dbcsr_scale(matrix_p, 0.5_dp) 397 398 DO istep = 1, admm_dm%mcweeny_max_steps 399 ! allocate new element in linked list 400 ALLOCATE (new_hist_entry) 401 new_hist_entry%next => history 402 history => new_hist_entry 403 history%count = istep 404 NULLIFY (new_hist_entry) 405 CALL dbcsr_create(history%m, template=matrix_p, matrix_type="N") 406 CALL dbcsr_copy(history%m, matrix_p, name="P from McWeeny") 407 408 ! calc PS and PSP 409 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, & 410 0.0_dp, matrix_ps) 411 412 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ps, matrix_p, & 413 0.0_dp, matrix_psp) 414 415 !test convergence 416 CALL dbcsr_copy(matrix_test, matrix_psp) 417 CALL dbcsr_add(matrix_test, matrix_p, 1.0_dp, -1.0_dp) 418 frob_norm = dbcsr_frobenius_norm(matrix_test) 419 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5,a,f16.8)') "McWeeny-Step", istep, & 420 ": Deviation of idempotency", frob_norm 421 IF (frob_norm < 1000_dp*admm_dm%eps_filter .AND. istep > 1) EXIT 422 423 ! build next P matrix 424 CALL dbcsr_copy(matrix_p, matrix_PSP, name="P from McWeeny") 425 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_PS, matrix_PSP, & 426 3.0_dp, matrix_p) 427 END DO 428 admm_dm%mcweeny_history(ispin)%p => history 429 IF (nspins == 1) CALL dbcsr_scale(matrix_p, 2.0_dp) 430 END DO 431 432 ! clean up 433 CALL dbcsr_release(matrix_PS) 434 CALL dbcsr_release(matrix_PSP) 435 CALL dbcsr_release(matrix_test) 436 CALL timestop(handle) 437 END SUBROUTINE purify_mcweeny 438 439! ************************************************************************************************** 440!> \brief Prepare auxilary KS-matrix for merge using reverse McWeeny 441!> \param ks_env ... 442!> \param matrix_ks_merge Output: The KS matrix for the merge 443!> \author Ole Schuett 444! ************************************************************************************************** 445 SUBROUTINE revert_purify_mcweeny(ks_env, matrix_ks_merge) 446 TYPE(qs_ks_env_type), POINTER :: ks_env 447 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge 448 449 CHARACTER(LEN=*), PARAMETER :: routineN = 'revert_purify_mcweeny', & 450 routineP = moduleN//':'//routineN 451 452 INTEGER :: handle, ispin, nspins, unit_nr 453 TYPE(admm_dm_type), POINTER :: admm_dm 454 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, & 455 matrix_s_aux_fit, & 456 matrix_s_aux_fit_vs_orb 457 TYPE(dbcsr_type), POINTER :: matrix_k 458 TYPE(dft_control_type), POINTER :: dft_control 459 TYPE(mcweeny_history_type), POINTER :: history_curr, history_next 460 461 CALL timeset(routineN, handle) 462 unit_nr = cp_logger_get_default_unit_nr() 463 NULLIFY (admm_dm, dft_control, matrix_ks, matrix_ks_aux_fit, & 464 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, & 465 history_next, history_curr, matrix_k) 466 467 CALL get_ks_env(ks_env, & 468 admm_dm=admm_dm, & 469 dft_control=dft_control, & 470 matrix_ks=matrix_ks, & 471 matrix_ks_aux_fit=matrix_ks_aux_fit, & 472 matrix_s_aux_fit=matrix_s_aux_fit, & 473 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb) 474 475 nspins = dft_control%nspins 476 ALLOCATE (matrix_ks_merge(nspins)) 477 478 DO ispin = 1, nspins 479 ALLOCATE (matrix_ks_merge(ispin)%matrix) 480 matrix_k => matrix_ks_merge(ispin)%matrix 481 CALL dbcsr_copy(matrix_k, matrix_ks_aux_fit(ispin)%matrix, name="K") 482 history_curr => admm_dm%mcweeny_history(ispin)%p 483 NULLIFY (admm_dm%mcweeny_history(ispin)%p) 484 485 ! reverse McWeeny iteration 486 DO WHILE (ASSOCIATED(history_curr)) 487 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5)') "Reverse McWeeny-Step ", history_curr%count 488 CALL reverse_mcweeny_step(matrix_k=matrix_k, & 489 matrix_s=matrix_s_aux_fit(1)%matrix, & 490 matrix_p=history_curr%m) 491 CALL dbcsr_release(history_curr%m) 492 history_next => history_curr%next 493 DEALLOCATE (history_curr) 494 history_curr => history_next 495 NULLIFY (history_next) 496 END DO 497 498 END DO 499 500 ! clean up 501 CALL timestop(handle) 502 503 END SUBROUTINE revert_purify_mcweeny 504 505! ************************************************************************************************** 506!> \brief Multiply matrix_k with partial derivative of McWeeny by reversing it. 507!> \param matrix_k ... 508!> \param matrix_s ... 509!> \param matrix_p ... 510!> \author Ole Schuett 511! ************************************************************************************************** 512 SUBROUTINE reverse_mcweeny_step(matrix_k, matrix_s, matrix_p) 513 TYPE(dbcsr_type) :: matrix_k, matrix_s, matrix_p 514 515 CHARACTER(LEN=*), PARAMETER :: routineN = 'reverse_mcweeny_step', & 516 routineP = moduleN//':'//routineN 517 518 INTEGER :: handle 519 TYPE(dbcsr_type) :: matrix_ps, matrix_sp, matrix_sum, & 520 matrix_tmp 521 522 CALL timeset(routineN, handle) 523 CALL dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N") 524 CALL dbcsr_create(matrix_sp, template=matrix_p, matrix_type="N") 525 CALL dbcsr_create(matrix_tmp, template=matrix_p, matrix_type="N") 526 CALL dbcsr_create(matrix_sum, template=matrix_p, matrix_type="N") 527 528 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, & 529 0.0_dp, matrix_ps) 530 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_p, & 531 0.0_dp, matrix_sp) 532 533 !TODO: can we exploid more symmetry? 534 CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_k, matrix_ps, & 535 0.0_dp, matrix_sum) 536 CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_sp, matrix_k, & 537 1.0_dp, matrix_sum) 538 539 !matrix_tmp = KPS 540 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_k, matrix_ps, & 541 0.0_dp, matrix_tmp) 542 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_tmp, matrix_ps, & 543 1.0_dp, matrix_sum) 544 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, & 545 1.0_dp, matrix_sum) 546 547 !matrix_tmp = SPK 548 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sp, matrix_k, & 549 0.0_dp, matrix_tmp) 550 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, & 551 1.0_dp, matrix_sum) 552 553 ! overwrite matrix_k 554 CALL dbcsr_copy(matrix_k, matrix_sum, name="K from reverse McWeeny") 555 556 ! clean up 557 CALL dbcsr_release(matrix_sum) 558 CALL dbcsr_release(matrix_tmp) 559 CALL dbcsr_release(matrix_ps) 560 CALL dbcsr_release(matrix_sp) 561 CALL timestop(handle) 562 END SUBROUTINE reverse_mcweeny_step 563 564END MODULE admm_dm_methods 565