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 to apply a delta pulse for RTP and EMD 8! ************************************************************************************************** 9 10MODULE rt_delta_pulse 11 USE cell_types, ONLY: cell_type 12 USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,& 13 cp_cfm_gemm 14 USE cp_cfm_diag, ONLY: cp_cfm_heevd 15 USE cp_cfm_types, ONLY: cp_cfm_create,& 16 cp_cfm_release,& 17 cp_cfm_to_cfm,& 18 cp_cfm_type 19 USE cp_control_types, ONLY: dft_control_type 20 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 21 copy_fm_to_dbcsr,& 22 cp_dbcsr_sm_fm_multiply 23 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,& 24 cp_fm_upper_to_full 25 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 26 cp_fm_cholesky_invert,& 27 cp_fm_cholesky_reduce,& 28 cp_fm_cholesky_restore 29 USE cp_fm_diag, ONLY: cp_fm_syevd 30 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 31 cp_fm_struct_release,& 32 cp_fm_struct_type 33 USE cp_fm_types, ONLY: cp_fm_create,& 34 cp_fm_get_info,& 35 cp_fm_p_type,& 36 cp_fm_release,& 37 cp_fm_set_all,& 38 cp_fm_to_fm,& 39 cp_fm_type 40 USE cp_gemm_interface, ONLY: cp_gemm 41 USE dbcsr_api, ONLY: dbcsr_copy,& 42 dbcsr_deallocate_matrix,& 43 dbcsr_filter,& 44 dbcsr_p_type,& 45 dbcsr_type 46 USE kinds, ONLY: dp 47 USE mathconstants, ONLY: twopi 48 USE qs_environment_types, ONLY: get_qs_env,& 49 qs_environment_type 50 USE qs_mo_types, ONLY: get_mo_set,& 51 mo_set_p_type 52 USE qs_moments, ONLY: build_berry_moment_matrix 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 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_delta_pulse' 62 63 PUBLIC :: apply_delta_pulse_periodic, & 64 apply_delta_pulse 65 66CONTAINS 67 68! ************************************************************************************************** 69!> \brief uses perturbation theory to get the proper initial conditions 70!> \param qs_env ... 71!> \param mos_old ... 72!> \param mos_new ... 73!> \author Joost & Martin (2011) 74! ************************************************************************************************** 75 76 SUBROUTINE apply_delta_pulse_periodic(qs_env, mos_old, mos_new) 77 TYPE(qs_environment_type), POINTER :: qs_env 78 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_old, mos_new 79 80 CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_periodic', & 81 routineP = moduleN//':'//routineN 82 83 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_sqrt 84 INTEGER :: handle, icol, idir, irow, ispin, nao, & 85 ncol_local, nmo, nrow_global, & 86 nrow_local, nvirt 87 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 88 REAL(KIND=dp) :: factor 89 REAL(KIND=dp), DIMENSION(3) :: kvec 90 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues 91 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_data 92 TYPE(cell_type), POINTER :: cell 93 TYPE(cp_cfm_type), POINTER :: oo_c, oo_v, oo_vt 94 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 95 TYPE(cp_fm_type), POINTER :: eigenvectors, mat_ks, mat_tmp, momentum, & 96 oo_1, oo_2, S_chol, S_inv_fm, tmpS, & 97 virtuals 98 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s 99 TYPE(dbcsr_type), POINTER :: S_inv 100 TYPE(dft_control_type), POINTER :: dft_control 101 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 102 TYPE(rt_prop_type), POINTER :: rtp 103 104 NULLIFY (dft_control) 105 CALL timeset(routineN, handle) 106 107 ! we need the overlap and ks matrix for a full diagionalization 108 CALL get_qs_env(qs_env, & 109 cell=cell, & 110 mos=mos, & 111 rtp=rtp, & 112 matrix_s=matrix_s, & 113 matrix_ks=matrix_ks, & 114 dft_control=dft_control) 115 CALL get_rtp(rtp=rtp, S_inv=S_inv) 116 CALL cp_fm_create(S_chol, matrix_struct=rtp%ao_ao_fmstruct, name="S_chol") 117 CALL cp_fm_create(S_inv_fm, matrix_struct=rtp%ao_ao_fmstruct, name="S_inv_fm") 118 CALL cp_fm_create(tmpS, matrix_struct=rtp%ao_ao_fmstruct) 119 CALL copy_dbcsr_to_fm(S_inv, S_inv_fm) 120 CALL cp_fm_upper_to_full(S_inv_fm, tmpS) 121 CALL cp_fm_get_info(S_inv_fm, nrow_global=nrow_global) 122 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol) 123 CALL cp_fm_cholesky_decompose(S_chol) 124 NULLIFY (mat_ks, eigenvectors, mat_tmp) 125 CALL cp_fm_create(mat_ks, matrix_struct=S_inv_fm%matrix_struct, name="mat_ks") 126 CALL cp_fm_create(eigenvectors, matrix_struct=S_inv_fm%matrix_struct, name="eigenvectors") 127 128 DO ispin = 1, SIZE(matrix_ks) 129 ALLOCATE (eigenvalues(nrow_global)) 130 CALL cp_fm_create(mat_tmp, matrix_struct=S_inv_fm%matrix_struct, name="mat_tmp") 131 132 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, mat_ks) 133 CALL cp_fm_cholesky_reduce(mat_ks, S_chol) 134 CALL cp_fm_syevd(mat_ks, mat_tmp, eigenvalues) 135 CALL cp_fm_cholesky_restore(mat_tmp, nrow_global, S_chol, eigenvectors, "SOLVE") 136 137 ! virtuals 138 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo) 139 nvirt = nao - nmo 140 CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, & 141 nrow_global=nrow_global, ncol_global=nvirt) 142 CALL cp_fm_create(virtuals, matrix_struct=fm_struct_tmp, name="virtuals") 143 CALL cp_fm_struct_release(fm_struct_tmp) 144 CALL cp_fm_to_fm(eigenvectors, virtuals, nvirt, nmo + 1, 1) 145 146 ! occupied 147 CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1)%matrix, nmo, 1, 1) 148 149 CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, & 150 nrow_global=nvirt, ncol_global=nmo) 151 CALL cp_fm_create(momentum, matrix_struct=fm_struct_tmp, name="momentum") 152 CALL cp_fm_struct_release(fm_struct_tmp) 153 154 ! the momentum operator (in a given direction) 155 CALL cp_fm_set_all(mos_new(2*ispin - 1)%matrix, 0.0_dp) 156 157 ! the prefactor (strength of the electric field) 158 kvec(:) = cell%h_inv(1, :)*dft_control%rtp_control%delta_pulse_direction(1) + & 159 cell%h_inv(2, :)*dft_control%rtp_control%delta_pulse_direction(2) + & 160 cell%h_inv(3, :)*dft_control%rtp_control%delta_pulse_direction(3) 161 kvec = -kvec*twopi*dft_control%rtp_control%delta_pulse_scale 162 163 DO idir = 1, 3 164 factor = kvec(idir) 165 IF (factor .NE. 0.0_dp) THEN 166 CALL cp_dbcsr_sm_fm_multiply(matrix_s(idir + 1)%matrix, mos_old(2*ispin - 1)%matrix, & 167 mos_old(2*ispin)%matrix, ncol=nmo) 168 CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1)%matrix, factor, mos_old(2*ispin)%matrix) 169 ENDIF 170 ENDDO 171 172 CALL cp_gemm('T', 'N', nvirt, nmo, nao, 1.0_dp, virtuals, mos_new(2*ispin - 1)%matrix, 0.0_dp, momentum) 173 174 ! the tricky bit ... rescale by the eigenvalue difference 175 CALL cp_fm_get_info(momentum, nrow_local=nrow_local, ncol_local=ncol_local, & 176 row_indices=row_indices, col_indices=col_indices, local_data=local_data) 177 DO icol = 1, ncol_local 178 DO irow = 1, nrow_local 179 factor = 1/(eigenvalues(col_indices(icol)) - eigenvalues(nmo + row_indices(irow))) 180 local_data(irow, icol) = factor*local_data(irow, icol) 181 ENDDO 182 ENDDO 183 CALL cp_fm_release(mat_tmp) 184 DEALLOCATE (eigenvalues) 185 186 ! now obtain the initial condition in mos_old 187 CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1)%matrix, nmo, 1, 1) 188 CALL cp_gemm("N", "N", nao, nmo, nvirt, 1.0_dp, virtuals, momentum, 0.0_dp, mos_old(2*ispin)%matrix) 189 190 CALL cp_fm_release(virtuals) 191 CALL cp_fm_release(momentum) 192 193 ! orthonormalize afterwards 194 CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, & 195 nrow_global=nmo, ncol_global=nmo) 196 CALL cp_fm_create(oo_1, matrix_struct=fm_struct_tmp, name="oo_1") 197 CALL cp_fm_create(oo_2, matrix_struct=fm_struct_tmp, name="oo_2") 198 CALL cp_fm_struct_release(fm_struct_tmp) 199 200 CALL cp_fm_create(mat_tmp, matrix_struct=mos_old(2*ispin - 1)%matrix%matrix_struct, name="tmp_mat") 201 ! get the complex overlap matrix 202 ! x^T S x + y^T S y + i (-y^TS x+x^T S y) 203 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mos_old(2*ispin - 1)%matrix, & 204 mat_tmp, ncol=nmo) 205 206 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*ispin - 1)%matrix, mat_tmp, 0.0_dp, oo_1) 207 CALL cp_gemm("T", "N", nmo, nmo, nao, -1.0_dp, mos_old(2*ispin)%matrix, mat_tmp, 0.0_dp, oo_2) 208 209 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mos_old(2*ispin)%matrix, & 210 mat_tmp, ncol=nmo) 211 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*ispin)%matrix, mat_tmp, 1.0_dp, oo_1) 212 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*ispin - 1)%matrix, mat_tmp, 1.0_dp, oo_2) 213 CALL cp_fm_release(mat_tmp) 214 215 CALL cp_cfm_create(oo_c, oo_1%matrix_struct) 216 CALL cp_cfm_create(oo_v, oo_1%matrix_struct) 217 CALL cp_cfm_create(oo_vt, oo_1%matrix_struct) 218 oo_c%local_data = CMPLX(oo_1%local_data, oo_2%local_data, KIND=dp) 219 220 ! compute inv(sqrt(overlap)) 221 ALLOCATE (eigenvalues(nmo)) 222 ALLOCATE (eigenvalues_sqrt(nmo)) 223 CALL cp_cfm_heevd(oo_c, oo_v, eigenvalues) 224 eigenvalues_sqrt = CMPLX(1.0_dp/SQRT(eigenvalues), 0.0_dp, dp) 225 CALL cp_cfm_to_cfm(oo_v, oo_vt) 226 CALL cp_cfm_column_scale(oo_v, eigenvalues_sqrt) 227 DEALLOCATE (eigenvalues) 228 DEALLOCATE (eigenvalues_sqrt) 229 CALL cp_cfm_gemm('N', 'C', nmo, nmo, nmo, (1.0_dp, 0.0_dp), & 230 oo_v, oo_vt, (0.0_dp, 0.0_dp), oo_c) 231 oo_1%local_data = REAL(oo_c%local_data, KIND=dp) 232 oo_2%local_data = AIMAG(oo_c%local_data) 233 CALL cp_cfm_release(oo_c) 234 CALL cp_cfm_release(oo_v) 235 CALL cp_cfm_release(oo_vt) 236 237 ! use this to compute the orthonormal vectors 238 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*ispin - 1)%matrix, oo_1, 0.0_dp, mos_new(2*ispin - 1)%matrix) 239 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*ispin - 1)%matrix, oo_2, 0.0_dp, mos_new(2*ispin)%matrix) 240 241 CALL cp_gemm("N", "N", nao, nmo, nmo, -1.0_dp, mos_old(2*ispin)%matrix, oo_2, 0.0_dp, mos_old(2*ispin - 1)%matrix) 242 CALL cp_fm_scale_and_add(1.0_dp, mos_old(2*ispin - 1)%matrix, 1.0_dp, mos_new(2*ispin - 1)%matrix) 243 244 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*ispin)%matrix, oo_1, 1.0_dp, mos_new(2*ispin)%matrix) 245 CALL cp_fm_to_fm(mos_new(2*ispin)%matrix, mos_old(2*ispin)%matrix) 246 247 CALL cp_fm_release(oo_1) 248 CALL cp_fm_release(oo_2) 249 END DO 250 251 CALL cp_fm_release(S_chol) 252 CALL cp_fm_release(mat_ks) 253 CALL cp_fm_release(eigenvectors) 254 255!*************************************************************** 256!remove later 257 CALL cp_fm_release(S_inv_fm) 258 CALL cp_fm_release(tmpS) 259!************************************************************** 260 CALL timestop(handle) 261 262 END SUBROUTINE apply_delta_pulse_periodic 263 264! ************************************************************************************************** 265!> \brief applies exp(ikr) to the wavefunction.... stored in mos_old... 266!> \param qs_env ... 267!> \param mos_old ... 268!> \param mos_new ... 269!> \author Joost & Martin (2011) 270! ************************************************************************************************** 271 272 SUBROUTINE apply_delta_pulse(qs_env, mos_old, mos_new) 273 TYPE(qs_environment_type), POINTER :: qs_env 274 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_old, mos_new 275 276 CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse', & 277 routineP = moduleN//':'//routineN 278 279 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_sqrt 280 INTEGER :: handle, i, nao, nmo 281 REAL(KIND=dp), DIMENSION(3) :: kvec 282 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues 283 TYPE(cell_type), POINTER :: cell 284 TYPE(cp_cfm_type), POINTER :: oo_c, oo_v, oo_vt 285 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 286 TYPE(cp_fm_type), POINTER :: mat_S, oo_1, oo_2, S_inv_fm, tmp 287 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 288 TYPE(dbcsr_type), POINTER :: cosmat, S_inv, sinmat 289 TYPE(dft_control_type), POINTER :: dft_control 290 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 291 TYPE(rt_prop_type), POINTER :: rtp 292 293 NULLIFY (dft_control) 294 295 CALL timeset(routineN, handle) 296 297 ! we need the inverse overlap 298 299 CALL get_qs_env(qs_env, & 300 mos=mos, & 301 rtp=rtp, & 302 matrix_s=matrix_s, & 303 dft_control=dft_control) 304 CALL get_rtp(rtp=rtp, S_inv=S_inv) 305 306 CALL cp_fm_create(S_inv_fm, matrix_struct=rtp%ao_ao_fmstruct, name="tmp_mat") 307 308 CALL cp_fm_create(tmp, matrix_struct=rtp%ao_ao_fmstruct, name="tmp_mat") 309 310 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_inv_fm) 311 CALL cp_fm_cholesky_decompose(S_inv_fm) 312 CALL cp_fm_cholesky_invert(S_inv_fm) 313 CALL cp_fm_upper_to_full(S_inv_fm, tmp) 314 315 CALL cp_fm_create(mat_S, matrix_struct=S_inv_fm%matrix_struct, name="mat_S") 316 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, mat_S) 317 CALL cp_fm_upper_to_full(mat_S, tmp) 318 319 CALL cp_fm_release(tmp) 320 321 ! we need the berry matrix 322 CALL get_qs_env(qs_env, cell=cell) 323 324 ! direction ... unscaled, this will yield a exp(ikr) that is periodic with the cell 325 kvec(:) = cell%h_inv(1, :)*dft_control%rtp_control%delta_pulse_direction(1) + & 326 cell%h_inv(2, :)*dft_control%rtp_control%delta_pulse_direction(2) + & 327 cell%h_inv(3, :)*dft_control%rtp_control%delta_pulse_direction(3) 328 kvec = -kvec*twopi 329 ! scaling will make the things not periodic with the cell, which would only be good for gas phase systems ? 330 kvec(:) = dft_control%rtp_control%delta_pulse_scale*kvec 331 332 ALLOCATE (cosmat, sinmat) 333 CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM') 334 CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM') 335 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec) 336 337 ! apply inv(S)*operator to C 338 DO i = 1, SIZE(mos) 339 CALL get_mo_set(mos(i)%mo_set, nao=nao, nmo=nmo) 340 CALL cp_dbcsr_sm_fm_multiply(cosmat, mos(i)%mo_set%mo_coeff, mos_new(2*i - 1)%matrix, ncol=nmo) 341 CALL cp_dbcsr_sm_fm_multiply(sinmat, mos(i)%mo_set%mo_coeff, mos_new(2*i)%matrix, ncol=nmo) 342 343 CALL cp_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i - 1)%matrix, 0.0_dp, mos_old(2*i - 1)%matrix) 344 CALL cp_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i)%matrix, 0.0_dp, mos_old(2*i)%matrix) 345 346 ! in a finite basis, unfortunately, inv(S)*operator is not unitary, so orthonormalize afterwards 347 CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, & 348 nrow_global=nmo, ncol_global=nmo) 349 CALL cp_fm_create(oo_1, matrix_struct=fm_struct_tmp, name="oo_1") 350 CALL cp_fm_create(oo_2, matrix_struct=fm_struct_tmp, name="oo_2") 351 CALL cp_fm_struct_release(fm_struct_tmp) 352 353 CALL cp_fm_create(tmp, matrix_struct=mos_old(2*i - 1)%matrix%matrix_struct, name="tmp_mat") 354 ! get the complex overlap matrix 355 ! x^T S x + y^T S y + i (-y^TS x+x^T S y) 356 CALL cp_gemm("N", "N", nao, nmo, nao, 1.0_dp, mat_S, mos_old(2*i - 1)%matrix, 0.0_dp, tmp) 357 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*i - 1)%matrix, tmp, 0.0_dp, oo_1) 358 CALL cp_gemm("T", "N", nmo, nmo, nao, -1.0_dp, mos_old(2*i)%matrix, tmp, 0.0_dp, oo_2) 359 360 CALL cp_gemm("N", "N", nao, nmo, nao, 1.0_dp, mat_S, mos_old(2*i)%matrix, 0.0_dp, tmp) 361 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*i)%matrix, tmp, 1.0_dp, oo_1) 362 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*i - 1)%matrix, tmp, 1.0_dp, oo_2) 363 CALL cp_fm_release(tmp) 364 365 CALL cp_cfm_create(oo_c, oo_1%matrix_struct) 366 CALL cp_cfm_create(oo_v, oo_1%matrix_struct) 367 CALL cp_cfm_create(oo_vt, oo_1%matrix_struct) 368 oo_c%local_data = CMPLX(oo_1%local_data, oo_2%local_data, KIND=dp) 369 370 ! compute inv(sqrt(overlap)) 371 ALLOCATE (eigenvalues(nmo)) 372 ALLOCATE (eigenvalues_sqrt(nmo)) 373 CALL cp_cfm_heevd(oo_c, oo_v, eigenvalues) 374 eigenvalues_sqrt = CMPLX(1.0_dp/SQRT(eigenvalues), 0.0_dp, dp) 375 CALL cp_cfm_to_cfm(oo_v, oo_vt) 376 CALL cp_cfm_column_scale(oo_v, eigenvalues_sqrt) 377 DEALLOCATE (eigenvalues) 378 DEALLOCATE (eigenvalues_sqrt) 379 CALL cp_cfm_gemm('N', 'C', nmo, nmo, nmo, (1.0_dp, 0.0_dp), & 380 oo_v, oo_vt, (0.0_dp, 0.0_dp), oo_c) 381 oo_1%local_data = REAL(oo_c%local_data, KIND=dp) 382 oo_2%local_data = AIMAG(oo_c%local_data) 383 CALL cp_cfm_release(oo_c) 384 CALL cp_cfm_release(oo_v) 385 CALL cp_cfm_release(oo_vt) 386 387 ! use this to compute the orthonormal vectors 388 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*i - 1)%matrix, oo_1, 0.0_dp, mos_new(2*i - 1)%matrix) 389 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*i - 1)%matrix, oo_2, 0.0_dp, mos_new(2*i)%matrix) 390 391 CALL cp_gemm("N", "N", nao, nmo, nmo, -1.0_dp, mos_old(2*i)%matrix, oo_2, 0.0_dp, mos_old(2*i - 1)%matrix) 392 CALL cp_fm_scale_and_add(1.0_dp, mos_old(2*i - 1)%matrix, 1.0_dp, mos_new(2*i - 1)%matrix) 393 394 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*i)%matrix, oo_1, 1.0_dp, mos_new(2*i)%matrix) 395 CALL cp_fm_to_fm(mos_new(2*i)%matrix, mos_old(2*i)%matrix) 396 397 CALL cp_fm_release(oo_1) 398 CALL cp_fm_release(oo_2) 399 END DO 400 401 CALL cp_fm_release(mat_S) 402 403 CALL dbcsr_deallocate_matrix(cosmat) 404 CALL dbcsr_deallocate_matrix(sinmat) 405 406!*************************************************************** 407!remove later 408 CALL copy_fm_to_dbcsr(S_inv_fm, S_inv) 409 CALL dbcsr_filter(S_inv, rtp%filter_eps) 410 CALL cp_fm_release(S_inv_fm) 411!************************************************************** 412 413 CALL timestop(handle) 414 415 END SUBROUTINE apply_delta_pulse 416 417END MODULE rt_delta_pulse 418