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