1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines needed for EMD 8!> \author Florian Schiffmann (02.09) 9! ************************************************************************************************** 10 11MODULE rt_propagation_forces 12 USE admm_types, ONLY: admm_type 13 USE atomic_kind_types, ONLY: atomic_kind_type,& 14 get_atomic_kind_set 15 USE cp_control_types, ONLY: dft_control_type,& 16 rtp_control_type 17 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 18 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,& 19 cp_dbcsr_sm_fm_multiply,& 20 dbcsr_create_dist_r_unrot 21 USE cp_fm_struct, ONLY: cp_fm_struct_type 22 USE cp_fm_types, ONLY: cp_fm_create,& 23 cp_fm_p_type 24 USE cp_fm_vect, ONLY: cp_fm_vect_dealloc 25 USE cp_gemm_interface, ONLY: cp_gemm 26 USE dbcsr_api, ONLY: & 27 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_release, & 28 dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, & 29 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 30 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry 31 USE kinds, ONLY: dp 32 USE mathconstants, ONLY: one,& 33 zero 34 USE particle_types, ONLY: particle_type 35 USE qs_environment_types, ONLY: get_qs_env,& 36 qs_environment_type 37 USE qs_force_types, ONLY: add_qs_force,& 38 qs_force_type 39 USE qs_ks_types, ONLY: qs_ks_env_type 40 USE qs_mo_types, ONLY: get_mo_set,& 41 mo_set_p_type 42 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 43 USE qs_overlap, ONLY: build_overlap_force 44 USE rt_propagation_types, ONLY: get_rtp,& 45 rt_prop_type 46#include "./base/base_uses.f90" 47 48 IMPLICIT NONE 49 PRIVATE 50 51 PUBLIC :: calc_c_mat_force, & 52 rt_admm_force 53 54 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_forces' 55 56CONTAINS 57 58! ************************************************************************************************** 59!> \brief calculates the three additional force contributions needed in EMD 60!> P_imag*C , P_imag*B*S^-1*S_der , P*S^-1*H*S_der 61!> driver routine 62!> \param qs_env ... 63!> \par History 64!> \author Florian Schiffmann (02.09) 65! ************************************************************************************************** 66 67 SUBROUTINE calc_c_mat_force(qs_env) 68 TYPE(qs_environment_type), POINTER :: qs_env 69 70 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_c_mat_force', & 71 routineP = moduleN//':'//routineN 72 73 IF (qs_env%rtp%linear_scaling) THEN 74 CALL calc_c_mat_force_ls(qs_env) 75 ELSE 76 CALL calc_c_mat_force_fm(qs_env) 77 END IF 78 79 END SUBROUTINE calc_c_mat_force 80 81! ************************************************************************************************** 82!> \brief standard treatment for fm MO based calculations 83!> P_imag*C , P_imag*B*S^-1*S_der , P*S^-1*H*S_der 84!> \param qs_env ... 85! ************************************************************************************************** 86 SUBROUTINE calc_c_mat_force_fm(qs_env) 87 TYPE(qs_environment_type), POINTER :: qs_env 88 89 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_c_mat_force_fm', & 90 routineP = moduleN//':'//routineN 91 92 INTEGER :: handle, i, im, ispin, nao, natom, nmo, re 93 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 94 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size 95 REAL(KIND=dp) :: alpha 96 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 97 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new 98 TYPE(dbcsr_distribution_type) :: dist, SinvB_dist 99 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat, S_der, SinvB, SinvH 100 TYPE(dbcsr_type) :: db_mo_tmp1, db_mo_tmp2, db_mos_im, & 101 db_mos_re 102 TYPE(dbcsr_type), POINTER :: rho_im_sparse, tmp_dbcsr 103 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 104 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 105 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 106 TYPE(rt_prop_type), POINTER :: rtp 107 108 CALL timeset(routineN, handle) 109 110 NULLIFY (rtp, particle_set, atomic_kind_set, mos) 111 NULLIFY (tmp_dbcsr, rho_im_sparse) 112 CALL get_qs_env(qs_env=qs_env, rtp=rtp, particle_set=particle_set, & 113 atomic_kind_set=atomic_kind_set, mos=mos, force=force) 114 115 CALL get_rtp(rtp=rtp, C_mat=C_mat, S_der=S_der, & 116 SinvH=SinvH, SinvB=SinvB, mos_new=mos_new) 117 118 ALLOCATE (tmp_dbcsr) 119 ALLOCATE (rho_im_sparse) 120 CALL dbcsr_create(tmp_dbcsr, template=SinvH(1)%matrix) 121 CALL dbcsr_create(rho_im_sparse, template=SinvH(1)%matrix) 122 123 natom = SIZE(particle_set) 124 ALLOCATE (atom_of_kind(natom)) 125 ALLOCATE (kind_of(natom)) 126 127 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of) 128 129 DO ispin = 1, SIZE(SinvH) 130 re = 2*ispin - 1 131 im = 2*ispin 132 alpha = mos(ispin)%mo_set%maxocc 133 134 CALL get_mo_set(mos(ispin)%mo_set, nao=nao, nmo=nmo) 135 136 NULLIFY (col_blk_size) 137 CALL dbcsr_get_info(SinvB(ispin)%matrix, distribution=SinvB_dist, row_blk_size=row_blk_size) 138 CALL dbcsr_create_dist_r_unrot(dist, SinvB_dist, nmo, col_blk_size) 139 140 CALL dbcsr_create(db_mos_re, "D", dist, dbcsr_type_no_symmetry, & 141 row_blk_size, col_blk_size, nze=0) 142 CALL dbcsr_create(db_mos_im, template=db_mos_re) 143 CALL dbcsr_create(db_mo_tmp1, template=db_mos_re) 144 CALL dbcsr_create(db_mo_tmp2, template=db_mos_re) 145 146 CALL copy_fm_to_dbcsr(mos_new(im)%matrix, db_mos_im) 147 CALL copy_fm_to_dbcsr(mos_new(re)%matrix, db_mos_re) 148 149 CALL dbcsr_multiply("N", "N", alpha, SinvB(ispin)%matrix, db_mos_im, 0.0_dp, db_mo_tmp1) 150 CALL dbcsr_multiply("N", "N", alpha, SinvH(ispin)%matrix, db_mos_re, 1.0_dp, db_mo_tmp1) 151 CALL dbcsr_multiply("N", "N", -alpha, SinvB(ispin)%matrix, db_mos_re, 0.0_dp, db_mo_tmp2) 152 CALL dbcsr_multiply("N", "N", alpha, SinvH(ispin)%matrix, db_mos_im, 1.0_dp, db_mo_tmp2) 153 CALL dbcsr_multiply("N", "T", 1.0_dp, db_mo_tmp1, db_mos_re, 0.0_dp, tmp_dbcsr) 154 CALL dbcsr_multiply("N", "T", 1.0_dp, db_mo_tmp2, db_mos_im, 1.0_dp, tmp_dbcsr) 155 156 CALL dbcsr_multiply("N", "T", alpha, db_mos_re, db_mos_im, 0.0_dp, rho_im_sparse) 157 CALL dbcsr_multiply("N", "T", -alpha, db_mos_im, db_mos_re, 1.0_dp, rho_im_sparse) 158 159 CALL compute_forces(force, tmp_dbcsr, S_der, rho_im_sparse, C_mat, kind_of, atom_of_kind) 160 161 CALL dbcsr_release(db_mos_re) 162 CALL dbcsr_release(db_mos_im) 163 CALL dbcsr_release(db_mo_tmp1) 164 CALL dbcsr_release(db_mo_tmp2) 165 166 DEALLOCATE (col_blk_size) 167 CALL dbcsr_distribution_release(dist) 168 169 END DO 170 171 DO i = 1, SIZE(force) 172 force(i)%ehrenfest(:, :) = -force(i)%ehrenfest(:, :) 173 END DO 174 175 CALL dbcsr_release(tmp_dbcsr) 176 CALL dbcsr_release(rho_im_sparse) 177 DEALLOCATE (atom_of_kind) 178 DEALLOCATE (kind_of) 179 DEALLOCATE (tmp_dbcsr) 180 DEALLOCATE (rho_im_sparse) 181 182 CALL timestop(handle) 183 184 END SUBROUTINE calc_c_mat_force_fm 185 186! ************************************************************************************************** 187!> \brief special treatment ofr linear scaling 188!> P_imag*C , P_imag*B*S^-1*S_der , P*S^-1*H*S_der 189!> \param qs_env ... 190!> \par History 191!> 02.2014 switched to dbcsr matrices [Samuel Andermatt] 192!> \author Florian Schiffmann (02.09) 193! ************************************************************************************************** 194 195 SUBROUTINE calc_c_mat_force_ls(qs_env) 196 TYPE(qs_environment_type), POINTER :: qs_env 197 198 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_c_mat_force_ls', & 199 routineP = moduleN//':'//routineN 200 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 201 202 INTEGER :: handle, i, im, ispin, natom, re 203 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 204 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 205 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat, rho_new, S_der, SinvB, SinvH 206 TYPE(dbcsr_type), POINTER :: S_inv, S_minus_half, tmp, tmp2, tmp3 207 TYPE(dft_control_type), POINTER :: dft_control 208 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 209 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 210 TYPE(rt_prop_type), POINTER :: rtp 211 TYPE(rtp_control_type), POINTER :: rtp_control 212 213 CALL timeset(routineN, handle) 214 NULLIFY (rtp, particle_set, atomic_kind_set, dft_control) 215 216 CALL get_qs_env(qs_env, & 217 rtp=rtp, & 218 particle_set=particle_set, & 219 atomic_kind_set=atomic_kind_set, & 220 force=force, & 221 dft_control=dft_control) 222 223 rtp_control => dft_control%rtp_control 224 CALL get_rtp(rtp=rtp, C_mat=C_mat, S_der=S_der, S_inv=S_inv, & 225 SinvH=SinvH, SinvB=SinvB) 226 227 natom = SIZE(particle_set) 228 ALLOCATE (atom_of_kind(natom)) 229 ALLOCATE (kind_of(natom)) 230 231 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of) 232 233 NULLIFY (tmp) 234 ALLOCATE (tmp) 235 CALL dbcsr_create(tmp, template=SinvB(1)%matrix) 236 NULLIFY (tmp2) 237 ALLOCATE (tmp2) 238 CALL dbcsr_create(tmp2, template=SinvB(1)%matrix) 239 NULLIFY (tmp3) 240 ALLOCATE (tmp3) 241 CALL dbcsr_create(tmp3, template=SinvB(1)%matrix) 242 243 CALL get_rtp(rtp=rtp, rho_new=rho_new, S_minus_half=S_minus_half) 244 245 DO ispin = 1, SIZE(SinvH) 246 re = 2*ispin - 1 247 im = 2*ispin 248 CALL dbcsr_multiply("N", "N", one, SinvH(ispin)%matrix, rho_new(re)%matrix, zero, tmp, & 249 filter_eps=rtp%filter_eps) 250 CALL dbcsr_multiply("N", "N", one, SinvB(ispin)%matrix, rho_new(im)%matrix, one, tmp, & 251 filter_eps=rtp%filter_eps) 252 253 CALL compute_forces(force, tmp, S_der, rho_new(im)%matrix, C_mat, kind_of, atom_of_kind) 254 255 END DO 256 257 ! recall QS forces, at this point have the other sign. 258 DO i = 1, SIZE(force) 259 force(i)%ehrenfest(:, :) = -force(i)%ehrenfest(:, :) 260 END DO 261 262 CALL dbcsr_deallocate_matrix(tmp) 263 CALL dbcsr_deallocate_matrix(tmp2) 264 CALL dbcsr_deallocate_matrix(tmp3) 265 266 DEALLOCATE (atom_of_kind) 267 DEALLOCATE (kind_of) 268 269 CALL timestop(handle) 270 271 END SUBROUTINE 272 273! ************************************************************************************************** 274!> \brief ... 275!> \param force ... 276!> \param tmp ... 277!> \param S_der ... 278!> \param rho_im ... 279!> \param C_mat ... 280!> \param kind_of ... 281!> \param atom_of_kind ... 282! ************************************************************************************************** 283 SUBROUTINE compute_forces(force, tmp, S_der, rho_im, C_mat, kind_of, atom_of_kind) 284 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 285 TYPE(dbcsr_type), POINTER :: tmp 286 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: S_der 287 TYPE(dbcsr_type), POINTER :: rho_im 288 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat 289 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, atom_of_kind 290 291 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_forces', routineP = moduleN//':'//routineN 292 293 INTEGER :: col_atom, i, ikind, kind_atom, row_atom 294 LOGICAL :: found 295 REAL(dp), DIMENSION(:), POINTER :: block_values, block_values2 296 TYPE(dbcsr_iterator_type) :: iter 297 298 DO i = 1, 3 299 !Calculate the sum over the hadmard product 300 !S_der part 301 302 CALL dbcsr_iterator_start(iter, tmp) 303 DO WHILE (dbcsr_iterator_blocks_left(iter)) 304 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values) 305 CALL dbcsr_get_block_p(S_der(i)%matrix, row_atom, col_atom, block_values2, found=found) 306 IF (found) THEN 307 ikind = kind_of(col_atom) 308 kind_atom = atom_of_kind(col_atom) 309 !The block_values are in a vector format, 310 ! so the dot_product is the sum over all elements of the hamand product, that I need 311 force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + & 312 2.0_dp*DOT_PRODUCT(block_values, block_values2) 313 ENDIF 314 END DO 315 CALL dbcsr_iterator_stop(iter) 316 317 !C_mat part 318 319 CALL dbcsr_iterator_start(iter, rho_im) 320 DO WHILE (dbcsr_iterator_blocks_left(iter)) 321 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values) 322 CALL dbcsr_get_block_p(C_mat(i)%matrix, row_atom, col_atom, block_values2, found=found) 323 IF (found) THEN 324 ikind = kind_of(col_atom) 325 kind_atom = atom_of_kind(col_atom) 326 !The block_values are in a vector format, so the dot_product is 327 ! the sum over all elements of the hamand product, that I need 328 force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + & 329 2.0_dp*DOT_PRODUCT(block_values, block_values2) 330 ENDIF 331 END DO 332 CALL dbcsr_iterator_stop(iter) 333 END DO 334 335 END SUBROUTINE compute_forces 336 337! ************************************************************************************************** 338!> \brief ... 339!> \param qs_env ... 340! ************************************************************************************************** 341 SUBROUTINE rt_admm_force(qs_env) 342 TYPE(qs_environment_type), POINTER :: qs_env 343 344 CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_admm_force', routineP = moduleN//':'//routineN 345 346 TYPE(admm_type), POINTER :: admm_env 347 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos, mos_admm 348 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_im, KS_aux_re, matrix_s_aux_fit, & 349 matrix_s_aux_fit_vs_orb 350 TYPE(rt_prop_type), POINTER :: rtp 351 352 CALL get_qs_env(qs_env, & 353 admm_env=admm_env, & 354 rtp=rtp, & 355 matrix_ks_aux_fit=KS_aux_re, & 356 matrix_ks_aux_fit_im=KS_aux_im, & 357 matrix_s_aux_fit=matrix_s_aux_fit, & 358 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb) 359 360 CALL get_rtp(rtp=rtp, mos_new=mos, admm_mos=mos_admm) 361 362 ! currently only none option 363 CALL rt_admm_forces_none(qs_env, admm_env, KS_aux_re, KS_aux_im, & 364 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_admm, mos) 365 366 END SUBROUTINE rt_admm_force 367 368! ************************************************************************************************** 369!> \brief ... 370!> \param qs_env ... 371!> \param admm_env ... 372!> \param KS_aux_re ... 373!> \param KS_aux_im ... 374!> \param matrix_s_aux_fit ... 375!> \param matrix_s_aux_fit_vs_orb ... 376!> \param mos_admm ... 377!> \param mos ... 378! ************************************************************************************************** 379 SUBROUTINE rt_admm_forces_none(qs_env, admm_env, KS_aux_re, KS_aux_im, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_admm, mos) 380 TYPE(qs_environment_type), POINTER :: qs_env 381 TYPE(admm_type), POINTER :: admm_env 382 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_re, KS_aux_im, matrix_s_aux_fit, & 383 matrix_s_aux_fit_vs_orb 384 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_admm, mos 385 386 CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_admm_forces_none', & 387 routineP = moduleN//':'//routineN 388 389 INTEGER :: im, ispin, nao, natom, naux, nmo, re 390 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force 391 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 392 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: tmp_aux_aux, tmp_aux_mo, tmp_aux_mo1, & 393 tmp_aux_nao 394 TYPE(cp_fm_struct_type), POINTER :: mstruct 395 TYPE(dbcsr_type), POINTER :: matrix_w_q, matrix_w_s 396 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 397 POINTER :: sab_aux_fit_asymm, sab_aux_fit_vs_orb 398 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 399 TYPE(qs_ks_env_type), POINTER :: ks_env 400 401 NULLIFY (sab_aux_fit_asymm, sab_aux_fit_vs_orb, ks_env) 402! CALL cp_fm_create(tmp_aux_aux,admm_env%fm_struct_tmp,name="fm matrix") 403 404 CALL get_qs_env(qs_env, & 405 sab_aux_fit_asymm=sab_aux_fit_asymm, & 406 sab_aux_fit_vs_orb=sab_aux_fit_vs_orb, & 407 ks_env=ks_env) 408 409 ALLOCATE (matrix_w_s) 410 CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, & 411 name='W MATRIX AUX S', matrix_type=dbcsr_type_no_symmetry) 412 CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, sab_aux_fit_asymm) 413 414 ALLOCATE (matrix_w_q) 415 CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, & 416 "W MATRIX AUX Q") 417 418 DO ispin = 1, SIZE(KS_aux_re) 419 re = 2*ispin - 1; im = 2*ispin 420 naux = admm_env%nao_aux_fit; nmo = admm_env%nmo(ispin); nao = admm_env%nao_orb 421 422 ALLOCATE (tmp_aux_aux(2), tmp_aux_nao(2), tmp_aux_mo(2), tmp_aux_mo1(2)) 423 CALL cp_fm_create(tmp_aux_aux(1)%matrix, admm_env%work_aux_aux%matrix_struct, name="taa") 424 CALL cp_fm_create(tmp_aux_aux(2)%matrix, admm_env%work_aux_aux%matrix_struct, name="taa") 425 CALL cp_fm_create(tmp_aux_nao(1)%matrix, admm_env%work_aux_orb%matrix_struct, name="tao") 426 CALL cp_fm_create(tmp_aux_nao(2)%matrix, admm_env%work_aux_orb%matrix_struct, name="tao") 427 mstruct => admm_env%work_aux_nmo(ispin)%matrix%matrix_struct 428 CALL cp_fm_create(tmp_aux_mo(1)%matrix, mstruct, name="tam") 429 CALL cp_fm_create(tmp_aux_mo(2)%matrix, mstruct, name="tam") 430 CALL cp_fm_create(tmp_aux_mo1(1)%matrix, mstruct, name="tam") 431 CALL cp_fm_create(tmp_aux_mo1(2)%matrix, mstruct, name="tam") 432 433! First calculate H=KS_aux*C~, real part ends on work_aux_aux2, imaginary part ends at work_aux_aux3 434 CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(re)%matrix, tmp_aux_mo(re)%matrix, nmo, 4.0_dp, 0.0_dp) 435 CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(im)%matrix, tmp_aux_mo(im)%matrix, nmo, 4.0_dp, 0.0_dp) 436 CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(im)%matrix, tmp_aux_mo(re)%matrix, nmo, -4.0_dp, 1.0_dp) 437 CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(re)%matrix, tmp_aux_mo(im)%matrix, nmo, 4.0_dp, 1.0_dp) 438 439! Next step compute S-1*H 440 CALL cp_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(re)%matrix, 0.0_dp, tmp_aux_mo1(re)%matrix) 441 CALL cp_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(im)%matrix, 0.0_dp, tmp_aux_mo1(im)%matrix) 442 443! Here we go on with Ws=S-1*H * C^H (take care of sign of the imaginary part!!!) 444 445 CALL cp_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(re)%matrix, mos(re)%matrix, 0.0_dp, & 446 tmp_aux_nao(re)%matrix) 447 CALL cp_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im)%matrix, mos(im)%matrix, 1.0_dp, & 448 tmp_aux_nao(re)%matrix) 449 CALL cp_gemm("N", "T", naux, nao, nmo, 1.0_dp, tmp_aux_mo1(re)%matrix, mos(im)%matrix, 0.0_dp, & 450 tmp_aux_nao(im)%matrix) 451 CALL cp_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im)%matrix, mos(re)%matrix, 1.0_dp, & 452 tmp_aux_nao(im)%matrix) 453 454! Let's do the final bit Wq=S-1*H * C^H * A^T 455 CALL cp_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(re)%matrix, admm_env%A, 0.0_dp, tmp_aux_aux(re)%matrix) 456 CALL cp_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(im)%matrix, admm_env%A, 0.0_dp, tmp_aux_aux(im)%matrix) 457 458 ! *** copy to sparse matrix 459 CALL copy_fm_to_dbcsr(tmp_aux_nao(re)%matrix, matrix_w_q, keep_sparsity=.TRUE.) 460 461 ! *** copy to sparse matrix 462 CALL copy_fm_to_dbcsr(tmp_aux_aux(re)%matrix, matrix_w_s, keep_sparsity=.TRUE.) 463 464! *** This can be done in one call w_total = w_alpha + w_beta 465 ! allocate force vector 466 CALL get_qs_env(qs_env=qs_env, natom=natom) 467 ALLOCATE (admm_force(3, natom)) 468 admm_force = 0.0_dp 469 CALL build_overlap_force(ks_env, admm_force, & 470 basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", & 471 sab_nl=sab_aux_fit_asymm, matrix_p=matrix_w_s) 472 CALL build_overlap_force(ks_env, admm_force, & 473 basis_type_a="AUX_FIT", basis_type_b="ORB", & 474 sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q) 475 ! add forces 476 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, & 477 force=force) 478 CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set) 479 DEALLOCATE (admm_force) 480 481 ! *** Deallocated weighted density matrices 482 CALL dbcsr_deallocate_matrix(matrix_w_s) 483 CALL dbcsr_deallocate_matrix(matrix_w_q) 484 CALL cp_fm_vect_dealloc(tmp_aux_aux) 485 CALL cp_fm_vect_dealloc(tmp_aux_nao) 486 CALL cp_fm_vect_dealloc(tmp_aux_mo) 487 CALL cp_fm_vect_dealloc(tmp_aux_mo1) 488 END DO 489 490 END SUBROUTINE rt_admm_forces_none 491 492END MODULE rt_propagation_forces 493