1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utilities for rtp in combination with admm methods 8!> adapted routines from admm_method (author Manuel Guidon) 9!> 10!> \par History Use new "force only" overlap routine [07.2014,JGH] 11!> \author Florian Schiffmann 12! ************************************************************************************************** 13MODULE rtp_admm_methods 14 USE admm_types, ONLY: admm_env_create,& 15 admm_type 16 USE cp_control_types, ONLY: admm_control_type,& 17 dft_control_type 18 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 19 copy_fm_to_dbcsr,& 20 cp_dbcsr_plus_fm_fm_t 21 USE cp_fm_basic_linalg, ONLY: cp_fm_upper_to_full 22 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 23 cp_fm_cholesky_invert 24 USE cp_fm_types, ONLY: cp_fm_get_info,& 25 cp_fm_p_type,& 26 cp_fm_to_fm,& 27 cp_fm_type 28 USE cp_gemm_interface, ONLY: cp_gemm 29 USE cp_para_types, ONLY: cp_para_env_type 30 USE dbcsr_api, ONLY: & 31 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, & 32 dbcsr_get_info, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry 33 USE hfx_admm_utils, ONLY: create_admm_xc_section 34 USE input_constants, ONLY: do_admm_basis_projection,& 35 do_admm_purify_none 36 USE input_section_types, ONLY: section_vals_get_subs_vals,& 37 section_vals_type 38 USE kinds, ONLY: dp 39 USE mathconstants, ONLY: one,& 40 zero 41 USE particle_types, ONLY: particle_type 42 USE pw_types, ONLY: pw_p_type 43 USE qs_collocate_density, ONLY: calculate_rho_elec 44 USE qs_environment_types, ONLY: get_qs_env,& 45 qs_environment_type,& 46 set_qs_env 47 USE qs_ks_types, ONLY: qs_ks_env_type 48 USE qs_mo_types, ONLY: get_mo_set,& 49 mo_set_p_type 50 USE qs_rho_types, ONLY: qs_rho_get,& 51 qs_rho_set,& 52 qs_rho_type 53 USE rt_propagation_types, ONLY: get_rtp,& 54 rt_prop_type 55 USE task_list_types, ONLY: task_list_type 56#include "./base/base_uses.f90" 57 58 IMPLICIT NONE 59 60 PRIVATE 61 62 ! *** Public subroutines *** 63 PUBLIC :: rtp_admm_calc_rho_aux, rtp_admm_merge_ks_matrix 64 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rtp_admm_methods' 66 67CONTAINS 68 69! ************************************************************************************************** 70!> \brief Compute the ADMM density matrix in case of rtp (complex MO's) 71!> 72!> \param qs_env ... 73!> \par History 74! ************************************************************************************************** 75 SUBROUTINE rtp_admm_calc_rho_aux(qs_env) 76 77 TYPE(qs_environment_type), POINTER :: qs_env 78 79 CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_calc_rho_aux', & 80 routineP = moduleN//':'//routineN 81 82 INTEGER :: handle, ispin, nspins 83 LOGICAL :: s_mstruct_changed 84 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_aux 85 TYPE(admm_type), POINTER :: admm_env 86 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit 87 TYPE(cp_para_env_type), POINTER :: para_env 88 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_aux, matrix_p_aux_im, & 89 matrix_s_aux_fit, & 90 matrix_s_aux_fit_vs_orb 91 TYPE(dft_control_type), POINTER :: dft_control 92 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit 93 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g_aux, rho_r_aux 94 TYPE(qs_ks_env_type), POINTER :: ks_env 95 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit 96 TYPE(rt_prop_type), POINTER :: rtp 97 TYPE(task_list_type), POINTER :: task_list_aux_fit 98 99 CALL timeset(routineN, handle) 100 NULLIFY (admm_env, matrix_p_aux, matrix_p_aux_im, mos, & 101 mos_aux_fit, para_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, & 102 ks_env, dft_control, tot_rho_r_aux, rho_r_aux, rho_g_aux, task_list_aux_fit) 103 104 CALL get_qs_env(qs_env, & 105 admm_env=admm_env, & 106 ks_env=ks_env, & 107 dft_control=dft_control, & 108 para_env=para_env, & 109 matrix_s_aux_fit=matrix_s_aux_fit, & 110 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, & 111 mos=mos, & 112 mos_aux_fit=mos_aux_fit, & 113 rtp=rtp, & 114 rho=rho, & 115 rho_aux_fit=rho_aux_fit, & 116 s_mstruct_changed=s_mstruct_changed, & 117 task_list_aux_fit=task_list_aux_fit) 118 119 IF (admm_env%do_gapw) THEN 120 CPABORT("GAPW ADMM not implemented for real time propagation") 121 END IF 122 123 nspins = dft_control%nspins 124 125 CALL get_rtp(rtp=rtp, admm_mos=rtp_coeff_aux_fit) 126 CALL rtp_admm_fit_mo_coeffs(qs_env, admm_env, dft_control%admm_control, para_env, & 127 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, & 128 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, & 129 s_mstruct_changed) 130 131 DO ispin = 1, nspins 132 CALL qs_rho_get(rho_aux_fit, & 133 rho_ao=matrix_p_aux, & 134 rho_ao_im=matrix_p_aux_im, & 135 rho_r=rho_r_aux, & 136 rho_g=rho_g_aux, & 137 tot_rho_r=tot_rho_r_aux) 138 139 CALL rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, & 140 matrix_p_aux(ispin)%matrix, & 141 matrix_p_aux_im(ispin)%matrix, & 142 ispin) 143 144 CALL calculate_rho_elec(matrix_p=matrix_p_aux(ispin)%matrix, & 145 rho=rho_r_aux(ispin), & 146 rho_gspace=rho_g_aux(ispin), & 147 total_rho=tot_rho_r_aux(ispin), & 148 ks_env=ks_env, soft_valid=.FALSE., & 149 basis_type="AUX_FIT", & 150 task_list_external=task_list_aux_fit) 151 END DO 152 CALL set_qs_env(qs_env, admm_env=admm_env) 153 CALL qs_rho_set(rho_aux_fit, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) 154 155 CALL timestop(handle) 156 157 END SUBROUTINE rtp_admm_calc_rho_aux 158 159! ************************************************************************************************** 160!> \brief ... 161!> \param admm_env ... 162!> \param rtp_coeff_aux_fit ... 163!> \param density_matrix_aux ... 164!> \param density_matrix_aux_im ... 165!> \param ispin ... 166! ************************************************************************************************** 167 SUBROUTINE rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, density_matrix_aux, & 168 density_matrix_aux_im, ispin) 169 TYPE(admm_type), POINTER :: admm_env 170 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit 171 TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im 172 INTEGER, INTENT(in) :: ispin 173 174 CHARACTER(len=*), PARAMETER :: routineN = 'rtp_admm_calculate_dm', & 175 routineP = moduleN//':'//routineN 176 177 INTEGER :: handle 178 179 CALL timeset(routineN, handle) 180 181 SELECT CASE (admm_env%purification_method) 182 CASE (do_admm_purify_none) 183 CALL calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, & 184 rtp_coeff_aux_fit, ispin) 185 CASE DEFAULT 186 CPWARN("only purification NONE possible with RTP/EMD at the moment") 187 END SELECT 188 189 CALL timestop(handle) 190 191 END SUBROUTINE rtp_admm_calculate_dm 192 193! ************************************************************************************************** 194!> \brief ... 195!> \param qs_env ... 196!> \param admm_env ... 197!> \param admm_control ... 198!> \param para_env ... 199!> \param matrix_s_aux_fit ... 200!> \param matrix_s_mixed ... 201!> \param mos ... 202!> \param mos_aux_fit ... 203!> \param rtp ... 204!> \param rtp_coeff_aux_fit ... 205!> \param geometry_did_change ... 206! ************************************************************************************************** 207 SUBROUTINE rtp_admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, & 208 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change) 209 210 TYPE(qs_environment_type), POINTER :: qs_env 211 TYPE(admm_type), POINTER :: admm_env 212 TYPE(admm_control_type), POINTER :: admm_control 213 TYPE(cp_para_env_type), POINTER :: para_env 214 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed 215 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit 216 TYPE(rt_prop_type), POINTER :: rtp 217 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit 218 LOGICAL, INTENT(IN) :: geometry_did_change 219 220 CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_fit_mo_coeffs', & 221 routineP = moduleN//':'//routineN 222 223 INTEGER :: handle, natoms 224 LOGICAL :: recalc_S 225 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 226 TYPE(section_vals_type), POINTER :: input, xc_section 227 228 CALL timeset(routineN, handle) 229 230 NULLIFY (xc_section, particle_set) 231 232 IF (.NOT. (ASSOCIATED(admm_env))) THEN 233 ! setup admm environment 234 CALL get_qs_env(qs_env, input=input, particle_set=particle_set) 235 natoms = SIZE(particle_set, 1) 236 CALL admm_env_create(admm_env, admm_control, mos, mos_aux_fit, & 237 para_env, natoms) 238 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 239 CALL create_admm_xc_section(qs_env=qs_env, xc_section=xc_section, & 240 admm_env=admm_env) 241 242 IF (admm_control%method /= do_admm_basis_projection) THEN 243 CPWARN("RTP requires BASIS_PROJECTION.") 244 END IF 245 END IF 246 247 recalc_S = geometry_did_change .OR. (rtp%iter == 0 .AND. (rtp%istep == rtp%i_start)) 248 249 SELECT CASE (admm_env%purification_method) 250 CASE (do_admm_purify_none) 251 CALL rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 252 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, recalc_S) 253 CASE DEFAULT 254 CPWARN("Purification method not implemented in combination with RTP") 255 END SELECT 256 257 CALL timestop(handle) 258 259 END SUBROUTINE rtp_admm_fit_mo_coeffs 260! ************************************************************************************************** 261!> \brief Calculates the MO coefficients for the auxiliary fitting basis set 262!> by minimizing int (psi_i - psi_aux_i)^2 using Lagrangian Multipliers 263!> 264!> \param qs_env ... 265!> \param admm_env The ADMM env 266!> \param para_env The parallel env 267!> \param matrix_s_aux_fit the overlap matrix of the auxiliary fitting basis set 268!> \param matrix_s_mixed the mixed overlap matrix of the auxiliary fitting basis 269!> set and the orbital basis set 270!> \param mos the MO's of the orbital basis set 271!> \param mos_aux_fit the MO's of the auxiliary fitting basis set 272!> \param rtp ... 273!> \param rtp_coeff_aux_fit ... 274!> \param geometry_did_change flag to indicate if the geomtry changed 275!> \par History 276!> 05.2008 created [Manuel Guidon] 277!> \author Manuel Guidon 278! ************************************************************************************************** 279 SUBROUTINE rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 280 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change) 281 282 TYPE(qs_environment_type), POINTER :: qs_env 283 TYPE(admm_type), POINTER :: admm_env 284 TYPE(cp_para_env_type), POINTER :: para_env 285 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed 286 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit 287 TYPE(rt_prop_type), POINTER :: rtp 288 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit 289 LOGICAL, INTENT(IN) :: geometry_did_change 290 291 CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_fit_mo_coeffs_none', & 292 routineP = moduleN//':'//routineN 293 294 INTEGER :: handle, ispin, nao_aux_fit, nao_orb, & 295 natoms, nmo, nmo_mos, nspins 296 REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux 297 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new 298 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 299 TYPE(dft_control_type), POINTER :: dft_control 300 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 301 TYPE(section_vals_type), POINTER :: input, xc_section 302 303 CALL timeset(routineN, handle) 304 305 NULLIFY (dft_control, particle_set) 306 307 IF (.NOT. (ASSOCIATED(admm_env))) THEN 308 CALL get_qs_env(qs_env, input=input, particle_set=particle_set, dft_control=dft_control) 309 natoms = SIZE(particle_set, 1) 310 CALL admm_env_create(admm_env, dft_control%admm_control, mos, mos_aux_fit, para_env, natoms) 311 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 312 CALL create_admm_xc_section(qs_env=qs_env, xc_section=xc_section, & 313 admm_env=admm_env) 314 END IF 315 316 nao_aux_fit = admm_env%nao_aux_fit 317 nao_orb = admm_env%nao_orb 318 nspins = SIZE(mos) 319 320 ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed 321 322 IF (geometry_did_change) THEN 323 CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv) 324 CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux) 325 CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S) 326 327 CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q) 328 329 !! Calculate S'_inverse 330 CALL cp_fm_cholesky_decompose(admm_env%S_inv) 331 CALL cp_fm_cholesky_invert(admm_env%S_inv) 332 !! Symmetrize the guy 333 CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux) 334 !! Calculate A=S'^(-1)*P 335 CALL cp_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, & 336 1.0_dp, admm_env%S_inv, admm_env%Q, 0.0_dp, & 337 admm_env%A) 338 END IF 339 340 ! *** Calculate the mo_coeffs for the fitting basis 341 DO ispin = 1, nspins 342 nmo = admm_env%nmo(ispin) 343 IF (nmo == 0) CYCLE 344 !! Lambda = C^(T)*B*C 345 CALL get_rtp(rtp=rtp, mos_new=mos_new) 346 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, occupation_numbers=occ_num, nmo=nmo_mos) 347 CALL get_mo_set(mos_aux_fit(ispin)%mo_set, mo_coeff=mo_coeff_aux_fit, & 348 occupation_numbers=occ_num_aux) 349 350 CALL cp_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, & 351 1.0_dp, admm_env%A, mos_new(2*ispin - 1)%matrix, 0.0_dp, & 352 rtp_coeff_aux_fit(2*ispin - 1)%matrix) 353 CALL cp_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, & 354 1.0_dp, admm_env%A, mos_new(2*ispin)%matrix, 0.0_dp, & 355 rtp_coeff_aux_fit(2*ispin)%matrix) 356 357 CALL cp_fm_to_fm(rtp_coeff_aux_fit(2*ispin - 1)%matrix, mo_coeff_aux_fit) 358 END DO 359 360 CALL timestop(handle) 361 362 END SUBROUTINE rtp_fit_mo_coeffs_none 363 364! ************************************************************************************************** 365!> \brief ... 366!> \param density_matrix_aux ... 367!> \param density_matrix_aux_im ... 368!> \param rtp_coeff_aux_fit ... 369!> \param ispin ... 370! ************************************************************************************************** 371 SUBROUTINE calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, & 372 rtp_coeff_aux_fit, ispin) 373 374 TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im 375 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit 376 INTEGER, INTENT(in) :: ispin 377 378 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rtp_admm_density', & 379 routineP = moduleN//':'//routineN 380 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 381 382 INTEGER :: handle, im, ncol, re 383 REAL(KIND=dp) :: alpha 384 385 CALL timeset(routineN, handle) 386 387 re = 2*ispin - 1; im = 2*ispin 388 alpha = 3*one - REAL(SIZE(rtp_coeff_aux_fit)/2, dp) 389 CALL dbcsr_set(density_matrix_aux, zero) 390 CALL cp_fm_get_info(rtp_coeff_aux_fit(re)%matrix, ncol_global=ncol) 391 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, & 392 matrix_v=rtp_coeff_aux_fit(re)%matrix, & 393 matrix_g=rtp_coeff_aux_fit(re)%matrix, & 394 ncol=ncol, & 395 alpha=alpha) 396 397 ! It is actually complex conjugate but i*i=-1 therefore it must be added 398 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, & 399 matrix_v=rtp_coeff_aux_fit(im)%matrix, & 400 matrix_g=rtp_coeff_aux_fit(im)%matrix, & 401 ncol=ncol, & 402 alpha=alpha) 403 404! compute the imaginary part of the dm 405 CALL dbcsr_set(density_matrix_aux_im, zero) 406 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, & 407 matrix_v=rtp_coeff_aux_fit(im)%matrix, & 408 matrix_g=rtp_coeff_aux_fit(re)%matrix, & 409 ncol=ncol, & 410 alpha=alpha) 411 alpha = -alpha 412 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, & 413 matrix_v=rtp_coeff_aux_fit(re)%matrix, & 414 matrix_g=rtp_coeff_aux_fit(im)%matrix, & 415 ncol=ncol, & 416 alpha=alpha) 417 418 CALL timestop(handle) 419 420 END SUBROUTINE calculate_rtp_admm_density 421 422! ************************************************************************************************** 423!> \brief ... 424!> \param qs_env ... 425! ************************************************************************************************** 426 SUBROUTINE rtp_admm_merge_ks_matrix(qs_env) 427 TYPE(qs_environment_type), POINTER :: qs_env 428 429 CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_merge_ks_matrix', & 430 routineP = moduleN//':'//routineN 431 432 INTEGER :: handle, ispin 433 TYPE(admm_type), POINTER :: admm_env 434 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, & 435 matrix_ks_aux_fit_im, matrix_ks_im 436 TYPE(dft_control_type), POINTER :: dft_control 437 438 NULLIFY (admm_env, dft_control, & 439 matrix_ks, matrix_ks_im, matrix_ks_aux_fit, matrix_ks_aux_fit_im) 440 CALL timeset(routineN, handle) 441 442 CALL get_qs_env(qs_env, & 443 admm_env=admm_env, & 444 dft_control=dft_control, & 445 matrix_ks=matrix_ks, & 446 matrix_ks_im=matrix_ks_im, & 447 matrix_ks_aux_fit=matrix_ks_aux_fit, & 448 matrix_ks_aux_fit_im=matrix_ks_aux_fit_im) 449 450 DO ispin = 1, dft_control%nspins 451 452 SELECT CASE (admm_env%purification_method) 453 CASE (do_admm_purify_none) 454 CALL rt_merge_ks_matrix_none(ispin, admm_env, & 455 matrix_ks, matrix_ks_aux_fit) 456 CALL rt_merge_ks_matrix_none(ispin, admm_env, & 457 matrix_ks_im, matrix_ks_aux_fit_im) 458 CASE DEFAULT 459 CPWARN("only purification NONE possible with RTP/EMD at the moment") 460 END SELECT 461 462 ENDDO !spin loop 463 CALL timestop(handle) 464 465 END SUBROUTINE rtp_admm_merge_ks_matrix 466 467! ************************************************************************************************** 468!> \brief ... 469!> \param ispin ... 470!> \param admm_env ... 471!> \param matrix_ks ... 472!> \param matrix_ks_aux_fit ... 473! ************************************************************************************************** 474 SUBROUTINE rt_merge_ks_matrix_none(ispin, admm_env, & 475 matrix_ks, matrix_ks_aux_fit) 476 INTEGER, INTENT(IN) :: ispin 477 TYPE(admm_type), POINTER :: admm_env 478 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit 479 480 CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_merge_ks_matrix_none', & 481 routineP = moduleN//':'//routineN 482 483 CHARACTER :: matrix_type_fit 484 INTEGER :: handle, nao_aux_fit, nao_orb, nmo 485 INTEGER, SAVE :: counter = 0 486 TYPE(dbcsr_type) :: matrix_ks_nosym 487 TYPE(dbcsr_type), POINTER :: matrix_k_tilde 488 489 CALL timeset(routineN, handle) 490 491 counter = counter + 1 492 nao_aux_fit = admm_env%nao_aux_fit 493 nao_orb = admm_env%nao_orb 494 nmo = admm_env%nmo(ispin) 495 CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks_aux_fit(ispin)%matrix, & 496 matrix_type=dbcsr_type_no_symmetry) 497 CALL dbcsr_set(matrix_ks_nosym, 0.0_dp) 498 CALL dbcsr_desymmetrize(matrix_ks_aux_fit(ispin)%matrix, matrix_ks_nosym) 499 500 CALL copy_dbcsr_to_fm(matrix_ks_nosym, admm_env%K(ispin)%matrix) 501 502 !! K*A 503 CALL cp_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, & 504 1.0_dp, admm_env%K(ispin)%matrix, admm_env%A, 0.0_dp, & 505 admm_env%work_aux_orb) 506 !! A^T*K*A 507 CALL cp_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, & 508 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, & 509 admm_env%work_orb_orb) 510 511 CALL dbcsr_get_info(matrix_ks_aux_fit(ispin)%matrix, matrix_type=matrix_type_fit) 512 513 NULLIFY (matrix_k_tilde) 514 ALLOCATE (matrix_k_tilde) 515 CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, & 516 name='MATRIX K_tilde', matrix_type=matrix_type_fit) 517 518 CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix) 519 CALL dbcsr_set(matrix_k_tilde, 0.0_dp) 520 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.) 521 522 CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp) 523 524 CALL dbcsr_deallocate_matrix(matrix_k_tilde) 525 CALL dbcsr_release(matrix_ks_nosym) 526 527 CALL timestop(handle) 528 529 END SUBROUTINE rt_merge_ks_matrix_none 530 531END MODULE rtp_admm_methods 532