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