1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> created 07.2005 9!> \author MI (07.2005) 10! ************************************************************************************************** 11MODULE qs_operators_ao 12 USE ai_moments, ONLY: diff_momop,& 13 diffop,& 14 moment 15 USE ai_os_rr, ONLY: os_rr_ovlp 16 USE basis_set_types, ONLY: gto_basis_set_p_type,& 17 gto_basis_set_type 18 USE block_p_types, ONLY: block_p_type 19 USE cell_types, ONLY: cell_type,& 20 pbc 21 USE cp_log_handling, ONLY: cp_get_default_logger,& 22 cp_logger_type 23 USE cp_para_types, ONLY: cp_para_env_type 24 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 25 dbcsr_get_matrix_type,& 26 dbcsr_p_type,& 27 dbcsr_type_antisymmetric,& 28 dbcsr_type_no_symmetry 29 USE kinds, ONLY: dp 30 USE mathconstants, ONLY: pi 31 USE orbital_pointers, ONLY: coset,& 32 init_orbital_pointers,& 33 ncoset 34 USE particle_types, ONLY: particle_type 35 USE qs_environment_types, ONLY: get_qs_env,& 36 qs_environment_type 37 USE qs_kind_types, ONLY: get_qs_kind,& 38 get_qs_kind_set,& 39 qs_kind_type 40 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 41 neighbor_list_iterate,& 42 neighbor_list_iterator_create,& 43 neighbor_list_iterator_p_type,& 44 neighbor_list_iterator_release,& 45 neighbor_list_set_p_type 46#include "./base/base_uses.f90" 47 48 IMPLICIT NONE 49 PRIVATE 50 51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_operators_ao' 52 53! *** Public subroutines *** 54 55 PUBLIC :: p_xyz_ao, rRc_xyz_ao, rRc_xyz_der_ao 56 PUBLIC :: build_lin_mom_matrix, build_ang_mom_matrix 57 58CONTAINS 59 60! ************************************************************************************************** 61!> \brief Calculation of the linear momentum matrix over 62!> Cartesian Gaussian functions. 63!> \param qs_env ... 64!> \param matrix ... 65!> \date 27.02.2009 66!> \author VW 67!> \version 1.0 68! ************************************************************************************************** 69 SUBROUTINE build_lin_mom_matrix(qs_env, matrix) 70 71 TYPE(qs_environment_type), POINTER :: qs_env 72 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix 73 74 CHARACTER(len=*), PARAMETER :: routineN = 'build_lin_mom_matrix', & 75 routineP = moduleN//':'//routineN 76 77 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, & 78 ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, nseta, nsetb, & 79 sgfa, sgfb 80 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 81 npgfb, nsgfa, nsgfb 82 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 83 LOGICAL :: found, new_atom_b 84 REAL(KIND=dp) :: dab, rab2 85 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work 86 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab, rr_work 87 REAL(KIND=dp), DIMENSION(3) :: rab 88 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 89 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 90 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral 91 TYPE(cell_type), POINTER :: cell 92 TYPE(cp_logger_type), POINTER :: logger 93 TYPE(cp_para_env_type), POINTER :: para_env 94 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 95 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 96 TYPE(neighbor_list_iterator_p_type), & 97 DIMENSION(:), POINTER :: nl_iterator 98 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 99 POINTER :: sab_orb 100 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 101 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 102 TYPE(qs_kind_type), POINTER :: qs_kind 103 104 CALL timeset(routineN, handle) 105 106 NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env) 107 NULLIFY (logger) 108 109 logger => cp_get_default_logger() 110 111 CALL get_qs_env(qs_env=qs_env, & 112 qs_kind_set=qs_kind_set, & 113 particle_set=particle_set, & 114 neighbor_list_id=neighbor_list_id, & 115 para_env=para_env, & 116 sab_orb=sab_orb, & 117 cell=cell) 118 119 nkind = SIZE(qs_kind_set) 120 natom = SIZE(particle_set) 121 122! *** Allocate work storage *** 123 124 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 125 maxco=maxco, & 126 maxlgto=maxlgto, & 127 maxsgf=maxsgf) 128 129 ldai = ncoset(maxlgto + 1) 130 CALL init_orbital_pointers(ldai) 131 132 ALLOCATE (rr_work(ldai, ldai, 3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3)) 133 rr_work(:, :, :) = 0.0_dp 134 intab(:, :, :) = 0.0_dp 135 work(:, :) = 0.0_dp 136 137 ALLOCATE (basis_set_list(nkind)) 138 DO ikind = 1, nkind 139 qs_kind => qs_kind_set(ikind) 140 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a) 141 IF (ASSOCIATED(basis_set_a)) THEN 142 basis_set_list(ikind)%gto_basis_set => basis_set_a 143 ELSE 144 NULLIFY (basis_set_list(ikind)%gto_basis_set) 145 END IF 146 END DO 147 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 148 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 149 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 150 iatom=iatom, jatom=jatom, r=rab) 151 basis_set_a => basis_set_list(ikind)%gto_basis_set 152 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 153 basis_set_b => basis_set_list(jkind)%gto_basis_set 154 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 155 ! basis ikind 156 first_sgfa => basis_set_a%first_sgf 157 la_max => basis_set_a%lmax 158 la_min => basis_set_a%lmin 159 npgfa => basis_set_a%npgf 160 nseta = basis_set_a%nset 161 nsgfa => basis_set_a%nsgf_set 162 rpgfa => basis_set_a%pgf_radius 163 set_radius_a => basis_set_a%set_radius 164 sphi_a => basis_set_a%sphi 165 zeta => basis_set_a%zet 166 ! basis jkind 167 first_sgfb => basis_set_b%first_sgf 168 lb_max => basis_set_b%lmax 169 lb_min => basis_set_b%lmin 170 npgfb => basis_set_b%npgf 171 nsetb = basis_set_b%nset 172 nsgfb => basis_set_b%nsgf_set 173 rpgfb => basis_set_b%pgf_radius 174 set_radius_b => basis_set_b%set_radius 175 sphi_b => basis_set_b%sphi 176 zetb => basis_set_b%zet 177 178 IF (inode == 1) last_jatom = 0 179 IF (jatom /= last_jatom) THEN 180 new_atom_b = .TRUE. 181 last_jatom = jatom 182 ELSE 183 new_atom_b = .FALSE. 184 END IF 185 186 IF (new_atom_b) THEN 187 IF (iatom <= jatom) THEN 188 irow = iatom 189 icol = jatom 190 ELSE 191 irow = jatom 192 icol = iatom 193 END IF 194 195 DO i = 1, 3 196 NULLIFY (integral(i)%block) 197 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, & 198 row=irow, col=icol, BLOCK=integral(i)%block, found=found) 199 CPASSERT(ASSOCIATED(INTEGRAL(i)%block)) 200 ENDDO 201 END IF 202 203 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 204 dab = SQRT(rab2) 205 206 DO iset = 1, nseta 207 208 ncoa = npgfa(iset)*ncoset(la_max(iset)) 209 sgfa = first_sgfa(1, iset) 210 211 DO jset = 1, nsetb 212 213 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 214 215 ncob = npgfb(jset)*ncoset(lb_max(jset)) 216 sgfb = first_sgfb(1, jset) 217 218 ! *** Calculate the primitive fermi contact integrals *** 219 220 CALL lin_mom(la_max(iset), la_min(iset), npgfa(iset), & 221 rpgfa(:, iset), zeta(:, iset), & 222 lb_max(jset), lb_min(jset), npgfb(jset), & 223 rpgfb(:, jset), zetb(:, jset), & 224 rab, intab, SIZE(rr_work, 1), rr_work) 225 226 ! *** Contraction step *** 227 228 DO i = 1, 3 229 230 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 231 1.0_dp, intab(1, 1, i), SIZE(intab, 1), & 232 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 233 0.0_dp, work(1, 1), SIZE(work, 1)) 234 235 IF (iatom <= jatom) THEN 236 237 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 238 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 239 work(1, 1), SIZE(work, 1), & 240 1.0_dp, integral(i)%block(sgfa, sgfb), & 241 SIZE(integral(i)%block, 1)) 242 243 ELSE 244 245 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, & 246 -1.0_dp, work(1, 1), SIZE(work, 1), & 247 sphi_a(1, sgfa), SIZE(sphi_a, 1), & 248 1.0_dp, integral(i)%block(sgfb, sgfa), & 249 SIZE(integral(i)%block, 1)) 250 251 ENDIF 252 253 ENDDO 254 255 ENDDO 256 257 ENDDO 258 259 END DO 260 CALL neighbor_list_iterator_release(nl_iterator) 261 262 ! *** Release work storage *** 263 264 DEALLOCATE (intab, work, integral, basis_set_list) 265 266! *** Print the spin orbit matrix, if requested *** 267 268 !IF (BTEST(cp_print_key_should_output(logger%iter_info,& 269 ! qs_env%input,"DFT%PRINT%AO_MATRICES/LINEAR_MOMENTUM"),cp_p_file)) THEN 270 ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/LINEA_MOMENTUM",& 271 ! extension=".Log") 272 ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw) 273 ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw) 274 ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw) 275 ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 276 ! "DFT%PRINT%AO_MATRICES/LINEAR_MOMENTUM") 277 !END IF 278 279 CALL timestop(handle) 280 281 END SUBROUTINE build_lin_mom_matrix 282 283! ************************************************************************************************** 284!> \brief Calculation of the primitive paramagnetic spin orbit integrals over 285!> Cartesian Gaussian-type functions. 286!> \param la_max ... 287!> \param la_min ... 288!> \param npgfa ... 289!> \param rpgfa ... 290!> \param zeta ... 291!> \param lb_max ... 292!> \param lb_min ... 293!> \param npgfb ... 294!> \param rpgfb ... 295!> \param zetb ... 296!> \param rab ... 297!> \param intab ... 298!> \param ldrr ... 299!> \param rr ... 300!> \date 02.03.2009 301!> \author VW 302!> \version 1.0 303! ************************************************************************************************** 304 SUBROUTINE lin_mom(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, & 305 rab, intab, ldrr, rr) 306 INTEGER, INTENT(IN) :: la_max, la_min, npgfa 307 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta 308 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb 309 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb 310 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 311 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: intab 312 INTEGER, INTENT(IN) :: ldrr 313 REAL(dp), DIMENSION(0:ldrr-1, 0:ldrr-1, 3), & 314 INTENT(INOUT) :: rr 315 316 CHARACTER(len=*), PARAMETER :: routineN = 'lin_mom', routineP = moduleN//':'//routineN 317 318 INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, & 319 ipgf, j, jpgf, la, lb, ma, mb, na, nb 320 REAL(dp) :: dab, dumx, dumy, dumz, f0, rab2, xhi, zet 321 REAL(dp), DIMENSION(3) :: rap, rbp 322 323! *** Calculate the distance of the centers a and c *** 324 325 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2 326 dab = SQRT(rab2) 327 328 ! *** Loop over all pairs of primitive Gaussian-type functions *** 329 330 na = 0 331 332 DO ipgf = 1, npgfa 333 334 nb = 0 335 336 DO jpgf = 1, npgfb 337 338 ! *** Screening *** 339 340 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN 341 DO j = nb + 1, nb + ncoset(lb_max) 342 DO i = na + 1, na + ncoset(la_max) 343 intab(i, j, 1) = 0.0_dp 344 intab(i, j, 2) = 0.0_dp 345 intab(i, j, 3) = 0.0_dp 346 ENDDO 347 ENDDO 348 nb = nb + ncoset(lb_max) 349 CYCLE 350 ENDIF 351 352 ! *** Calculate some prefactors *** 353 zet = zeta(ipgf) + zetb(jpgf) 354 xhi = zeta(ipgf)*zetb(jpgf)/zet 355 rap = zetb(jpgf)*rab/zet 356 rbp = -zeta(ipgf)*rab/zet 357 358 f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2) 359 360 ! *** Calculate the recurrence relation *** 361 362 CALL os_rr_ovlp(rap, la_max + 1, rbp, lb_max, zet, ldrr, rr) 363 364 ! *** Calculate the primitive linear momentum integrals *** 365 DO lb = lb_min, lb_max 366 DO bx = 0, lb 367 DO by = 0, lb - bx 368 bz = lb - bx - by 369 cob = coset(bx, by, bz) 370 mb = nb + cob 371 DO la = la_min, la_max 372 DO ax = 0, la 373 DO ay = 0, la - ax 374 az = la - ax - ay 375 coa = coset(ax, ay, az) 376 ma = na + coa 377 ! 378 ! 379 ! (a|p_x|b) = 2*a*(a+1x|b) - N_x(a)*(a-1_x|b) 380 dumx = 2.0_dp*zeta(ipgf)*rr(ax + 1, bx, 1) 381 IF (ax .GT. 0) dumx = dumx - REAL(ax, dp)*rr(ax - 1, bx, 1) 382 intab(ma, mb, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3) 383 ! 384 ! (a|p_y|b) 385 dumy = 2.0_dp*zeta(ipgf)*rr(ay + 1, by, 2) 386 IF (ay .GT. 0) dumy = dumy - REAL(ay, dp)*rr(ay - 1, by, 2) 387 intab(ma, mb, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3) 388 ! 389 ! (a|p_z|b) 390 dumz = 2.0_dp*zeta(ipgf)*rr(az + 1, bz, 3) 391 IF (az .GT. 0) dumz = dumz - REAL(az, dp)*rr(az - 1, bz, 3) 392 intab(ma, mb, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz 393 ! 394 ENDDO 395 ENDDO 396 ENDDO !la 397 398 ENDDO 399 ENDDO 400 ENDDO !lb 401 402 nb = nb + ncoset(lb_max) 403 404 ENDDO 405 406 na = na + ncoset(la_max) 407 408 ENDDO 409 410 END SUBROUTINE lin_mom 411 412! ************************************************************************************************** 413!> \brief Calculation of the angular momentum matrix over 414!> Cartesian Gaussian functions. 415!> \param qs_env ... 416!> \param matrix ... 417!> \param rc ... 418!> \date 27.02.2009 419!> \author VW 420!> \version 1.0 421! ************************************************************************************************** 422 423 SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc) 424 425 TYPE(qs_environment_type), POINTER :: qs_env 426 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix 427 REAL(dp), DIMENSION(:), INTENT(IN) :: rc 428 429 CHARACTER(len=*), PARAMETER :: routineN = 'build_ang_mom_matrix', & 430 routineP = moduleN//':'//routineN 431 432 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, & 433 ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, nseta, nsetb, & 434 sgfa, sgfb 435 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 436 npgfb, nsgfa, nsgfb 437 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 438 LOGICAL :: found, new_atom_b 439 REAL(KIND=dp) :: dab, rab2 440 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work 441 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab, rr_work 442 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rbc 443 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 444 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 445 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral 446 TYPE(cell_type), POINTER :: cell 447 TYPE(cp_logger_type), POINTER :: logger 448 TYPE(cp_para_env_type), POINTER :: para_env 449 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 450 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 451 TYPE(neighbor_list_iterator_p_type), & 452 DIMENSION(:), POINTER :: nl_iterator 453 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 454 POINTER :: sab_all 455 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 456 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 457 TYPE(qs_kind_type), POINTER :: qs_kind 458 459 CALL timeset(routineN, handle) 460 461 NULLIFY (cell, sab_all, qs_kind_set, particle_set, para_env) 462 NULLIFY (logger) 463 464 logger => cp_get_default_logger() 465 466 CALL get_qs_env(qs_env=qs_env, & 467 qs_kind_set=qs_kind_set, & 468 particle_set=particle_set, & 469 neighbor_list_id=neighbor_list_id, & 470 para_env=para_env, & 471 sab_all=sab_all, & 472 cell=cell) 473 474 nkind = SIZE(qs_kind_set) 475 natom = SIZE(particle_set) 476 477! *** Allocate work storage *** 478 479 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 480 maxco=maxco, & 481 maxlgto=maxlgto, & 482 maxsgf=maxsgf) 483 484 ldai = ncoset(maxlgto + 1) 485 CALL init_orbital_pointers(ldai) 486 487 ALLOCATE (rr_work(ldai, ldai, 3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3)) 488 rr_work(:, :, :) = 0.0_dp 489 intab(:, :, :) = 0.0_dp 490 work(:, :) = 0.0_dp 491 492 ALLOCATE (basis_set_list(nkind)) 493 DO ikind = 1, nkind 494 qs_kind => qs_kind_set(ikind) 495 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a) 496 IF (ASSOCIATED(basis_set_a)) THEN 497 basis_set_list(ikind)%gto_basis_set => basis_set_a 498 ELSE 499 NULLIFY (basis_set_list(ikind)%gto_basis_set) 500 END IF 501 END DO 502 CALL neighbor_list_iterator_create(nl_iterator, sab_all) 503 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 504 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 505 iatom=iatom, jatom=jatom, r=rab) 506 basis_set_a => basis_set_list(ikind)%gto_basis_set 507 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 508 basis_set_b => basis_set_list(jkind)%gto_basis_set 509 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 510 ra = pbc(particle_set(iatom)%r, cell) 511 ! basis ikind 512 first_sgfa => basis_set_a%first_sgf 513 la_max => basis_set_a%lmax 514 la_min => basis_set_a%lmin 515 npgfa => basis_set_a%npgf 516 nseta = basis_set_a%nset 517 nsgfa => basis_set_a%nsgf_set 518 rpgfa => basis_set_a%pgf_radius 519 set_radius_a => basis_set_a%set_radius 520 sphi_a => basis_set_a%sphi 521 zeta => basis_set_a%zet 522 ! basis jkind 523 first_sgfb => basis_set_b%first_sgf 524 lb_max => basis_set_b%lmax 525 lb_min => basis_set_b%lmin 526 npgfb => basis_set_b%npgf 527 nsetb = basis_set_b%nset 528 nsgfb => basis_set_b%nsgf_set 529 rpgfb => basis_set_b%pgf_radius 530 set_radius_b => basis_set_b%set_radius 531 sphi_b => basis_set_b%sphi 532 zetb => basis_set_b%zet 533 534 IF (inode == 1) last_jatom = 0 535 536 IF (jatom /= last_jatom) THEN 537 new_atom_b = .TRUE. 538 last_jatom = jatom 539 ELSE 540 new_atom_b = .FALSE. 541 END IF 542 543 IF (new_atom_b) THEN 544 !IF (iatom <= jatom) THEN 545 irow = iatom 546 icol = jatom 547 !ELSE 548 ! irow = jatom 549 ! icol = iatom 550 !END IF 551 552 DO i = 1, 3 553 NULLIFY (INTEGRAL(i)%block) 554 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, & 555 row=irow, col=icol, BLOCK=INTEGRAL(i)%block, found=found) 556 CPASSERT(ASSOCIATED(INTEGRAL(i)%block)) 557 ENDDO 558 END IF 559 560 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 561 dab = SQRT(rab2) 562 563 DO iset = 1, nseta 564 565 ncoa = npgfa(iset)*ncoset(la_max(iset)) 566 sgfa = first_sgfa(1, iset) 567 568 DO jset = 1, nsetb 569 570 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 571 572 !IF(PRESENT(wancen)) THEN 573 ! rc = wancen 574 rac = pbc(rc, ra, cell) 575 rbc = rac + rab 576 !ELSE 577 ! rc(1:3) = rb(1:3) 578 ! rac(1:3) = -rab(1:3) 579 ! rbc(1:3) = 0.0_dp 580 !ENDIF 581 582 ncob = npgfb(jset)*ncoset(lb_max(jset)) 583 sgfb = first_sgfb(1, jset) 584 585 ! *** Calculate the primitive angular momentum integrals *** 586 587 CALL ang_mom(la_max(iset), la_min(iset), npgfa(iset), & 588 rpgfa(:, iset), zeta(:, iset), & 589 lb_max(jset), lb_min(jset), npgfb(jset), & 590 rpgfb(:, jset), zetb(:, jset), & 591 rab, rac, intab, SIZE(rr_work, 1), rr_work) 592 593 ! *** Contraction step *** 594 595 DO i = 1, 3 596 597 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 598 1.0_dp, intab(1, 1, i), SIZE(intab, 1), & 599 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 600 0.0_dp, work(1, 1), SIZE(work, 1)) 601 602 !IF (iatom <= jatom) THEN 603 604 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 605 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 606 work(1, 1), SIZE(work, 1), & 607 1.0_dp, integral(i)%block(sgfa, sgfb), & 608 SIZE(integral(i)%block, 1)) 609 610 !ELSE 611 ! 612 ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,& 613 ! -1.0_dp,work(1,1),SIZE(work,1),& 614 ! sphi_a(1,sgfa),SIZE(sphi_a,1),& 615 ! 1.0_dp,integral(i)%block(sgfb,sgfa),& 616 ! SIZE(integral(i)%block,1)) 617 ! 618 !ENDIF 619 620 ENDDO 621 622 ENDDO 623 624 ENDDO 625 626 END DO 627 CALL neighbor_list_iterator_release(nl_iterator) 628 629 ! *** Release work storage *** 630 631 DEALLOCATE (intab, work, integral, basis_set_list) 632 633! *** Print the spin orbit matrix, if requested *** 634 635 !IF (BTEST(cp_print_key_should_output(logger%iter_info,& 636 ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN 637 ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",& 638 ! extension=".Log") 639 ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw) 640 ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw) 641 ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw) 642 ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 643 ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM") 644 !END IF 645 646 CALL timestop(handle) 647 648 END SUBROUTINE build_ang_mom_matrix 649 650! ************************************************************************************************** 651!> \brief Calculation of the primitive paramagnetic spin orbit integrals over 652!> Cartesian Gaussian-type functions. 653!> \param la_max ... 654!> \param la_min ... 655!> \param npgfa ... 656!> \param rpgfa ... 657!> \param zeta ... 658!> \param lb_max ... 659!> \param lb_min ... 660!> \param npgfb ... 661!> \param rpgfb ... 662!> \param zetb ... 663!> \param rab ... 664!> \param rac ... 665!> \param intab ... 666!> \param ldrr ... 667!> \param rr ... 668!> \date 02.03.2009 669!> \author VW 670!> \version 1.0 671! ************************************************************************************************** 672 SUBROUTINE ang_mom(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, & 673 rab, rac, intab, ldrr, rr) 674 INTEGER, INTENT(IN) :: la_max, la_min, npgfa 675 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta 676 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb 677 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb 678 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab, rac 679 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: intab 680 INTEGER, INTENT(IN) :: ldrr 681 REAL(dp), DIMENSION(0:ldrr-1, 0:ldrr-1, 3), & 682 INTENT(INOUT) :: rr 683 684 CHARACTER(len=*), PARAMETER :: routineN = 'ang_mom', routineP = moduleN//':'//routineN 685 686 INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, & 687 ipgf, j, jpgf, la, lb, ma, mb, na, nb 688 REAL(dp) :: dab, dumx, dumy, dumz, f0, rab2, xhi, zet 689 REAL(dp), DIMENSION(3) :: rap, rbp 690 691! *** Calculate the distance of the centers a and c *** 692 693 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2 694 dab = SQRT(rab2) 695 696 ! *** Loop over all pairs of primitive Gaussian-type functions *** 697 698 na = 0 699 700 DO ipgf = 1, npgfa 701 702 nb = 0 703 704 DO jpgf = 1, npgfb 705 706 ! *** Screening *** 707 708 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN 709 DO j = nb + 1, nb + ncoset(lb_max) 710 DO i = na + 1, na + ncoset(la_max) 711 intab(i, j, 1) = 0.0_dp 712 intab(i, j, 2) = 0.0_dp 713 intab(i, j, 3) = 0.0_dp 714 ENDDO 715 ENDDO 716 nb = nb + ncoset(lb_max) 717 CYCLE 718 ENDIF 719 720 ! *** Calculate some prefactors *** 721 zet = zeta(ipgf) + zetb(jpgf) 722 xhi = zeta(ipgf)*zetb(jpgf)/zet 723 rap = zetb(jpgf)*rab/zet 724 rbp = -zeta(ipgf)*rab/zet 725 726 f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2) 727 728 ! *** Calculate the recurrence relation *** 729 730 CALL os_rr_ovlp(rap, la_max + 1, rbp, lb_max, zet, ldrr, rr) 731 732 ! *** Calculate the primitive Fermi contact integrals *** 733 734 DO lb = lb_min, lb_max 735 DO bx = 0, lb 736 DO by = 0, lb - bx 737 bz = lb - bx - by 738 cob = coset(bx, by, bz) 739 mb = nb + cob 740 DO la = la_min, la_max 741 DO ax = 0, la 742 DO ay = 0, la - ax 743 az = la - ax - ay 744 coa = coset(ax, ay, az) 745 ma = na + coa 746 ! 747 dumx = -2.0_dp*zeta(ipgf)*rr(ax + 1, bx, 1) 748 dumy = -2.0_dp*zeta(ipgf)*rr(ay + 1, by, 2) 749 dumz = -2.0_dp*zeta(ipgf)*rr(az + 1, bz, 3) 750 IF (ax .GT. 0) dumx = dumx + REAL(ax, dp)*rr(ax - 1, bx, 1) 751 IF (ay .GT. 0) dumy = dumy + REAL(ay, dp)*rr(ay - 1, by, 2) 752 IF (az .GT. 0) dumz = dumz + REAL(az, dp)*rr(az - 1, bz, 3) 753 ! 754 ! (a|l_z|b) 755 intab(ma, mb, 1) = -f0*rr(ax, bx, 1)*( & 756 & (rr(ay + 1, by, 2) + rac(2)*rr(ay, by, 2))*dumz & 757 & - (rr(az + 1, bz, 3) + rac(3)*rr(az, bz, 3))*dumy) 758 ! 759 ! (a|l_y|b) 760 intab(ma, mb, 2) = -f0*rr(ay, by, 2)*( & 761 & (rr(az + 1, bz, 3) + rac(3)*rr(az, bz, 3))*dumx & 762 & - (rr(ax + 1, bx, 1) + rac(1)*rr(ax, bx, 1))*dumz) 763 ! 764 ! (a|l_z|b) 765 intab(ma, mb, 3) = -f0*rr(az, bz, 3)*( & 766 & (rr(ax + 1, bx, 1) + rac(1)*rr(ax, bx, 1))*dumy & 767 & - (rr(ay + 1, by, 2) + rac(2)*rr(ay, by, 2))*dumx) 768 ! 769 ENDDO 770 ENDDO 771 ENDDO !la 772 773 ENDDO 774 ENDDO 775 ENDDO !lb 776 777 nb = nb + ncoset(lb_max) 778 779 ENDDO 780 781 na = na + ncoset(la_max) 782 783 ENDDO 784 785 END SUBROUTINE ang_mom 786 787! ************************************************************************************************** 788!> \brief Calculation of the components of the dipole operator in the velocity form 789!> The elements of the sparse matrices are the integrals in the 790!> basis functions 791!> \param op matrix representation of the p operator 792!> calculated in terms of the contracted basis functions 793!> \param qs_env environment for the lists and the basis sets 794!> \param minimum_image take into account only the first neighbors in the lists 795!> \par History 796!> 06.2005 created [MI] 797!> \author MI 798! ************************************************************************************************** 799 800 SUBROUTINE p_xyz_ao(op, qs_env, minimum_image) 801 802 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op 803 TYPE(qs_environment_type), POINTER :: qs_env 804 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image 805 806 CHARACTER(len=*), PARAMETER :: routineN = 'p_xyz_ao', routineP = moduleN//':'//routineN 807 808 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, & 809 ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb 810 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 811 npgfb, nsgfa, nsgfb 812 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 813 LOGICAL :: found, my_minimum_image, new_atom_b 814 REAL(KIND=dp) :: alpha, dab, Lxo2, Lyo2, Lzo2, rab2 815 REAL(KIND=dp), DIMENSION(3) :: ra, rab 816 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 817 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, & 818 zeta, zetb 819 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: difab 820 TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip 821 TYPE(cell_type), POINTER :: cell 822 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 823 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 824 TYPE(neighbor_list_iterator_p_type), & 825 DIMENSION(:), POINTER :: nl_iterator 826 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 827 POINTER :: sab_orb 828 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 829 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 830 TYPE(qs_kind_type), POINTER :: qs_kind 831 832 CALL timeset(routineN, handle) 833 834 NULLIFY (qs_kind, qs_kind_set) 835 NULLIFY (cell, particle_set) 836 NULLIFY (sab_orb) 837 NULLIFY (difab, op_dip, work) 838 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb) 839 NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb) 840 841 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 842 cell=cell, particle_set=particle_set, & 843 sab_orb=sab_orb) 844 845 nkind = SIZE(qs_kind_set) 846 847 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 848 maxco=ldwork, maxlgto=maxl) 849 850 my_minimum_image = .FALSE. 851 IF (PRESENT(minimum_image)) THEN 852 my_minimum_image = minimum_image 853 Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp 854 Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp 855 Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp 856 END IF 857 858 ldab = ldwork 859 860 ALLOCATE (difab(ldab, ldab, 3)) 861 difab(1:ldab, 1:ldab, 1:3) = 0.0_dp 862 ALLOCATE (work(ldwork, ldwork)) 863 work(1:ldwork, 1:ldwork) = 0.0_dp 864 ALLOCATE (op_dip(3)) 865 866 DO i = 1, 3 867 NULLIFY (op_dip(i)%block) 868 END DO 869 870 ALLOCATE (basis_set_list(nkind)) 871 DO ikind = 1, nkind 872 qs_kind => qs_kind_set(ikind) 873 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a) 874 IF (ASSOCIATED(basis_set_a)) THEN 875 basis_set_list(ikind)%gto_basis_set => basis_set_a 876 ELSE 877 NULLIFY (basis_set_list(ikind)%gto_basis_set) 878 END IF 879 END DO 880 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 881 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 882 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 883 iatom=iatom, jatom=jatom, r=rab) 884 basis_set_a => basis_set_list(ikind)%gto_basis_set 885 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 886 basis_set_b => basis_set_list(jkind)%gto_basis_set 887 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 888 ra = pbc(particle_set(iatom)%r, cell) 889 ! basis ikind 890 first_sgfa => basis_set_a%first_sgf 891 la_max => basis_set_a%lmax 892 la_min => basis_set_a%lmin 893 npgfa => basis_set_a%npgf 894 nseta = basis_set_a%nset 895 nsgfa => basis_set_a%nsgf_set 896 rpgfa => basis_set_a%pgf_radius 897 set_radius_a => basis_set_a%set_radius 898 sphi_a => basis_set_a%sphi 899 zeta => basis_set_a%zet 900 ! basis jkind 901 first_sgfb => basis_set_b%first_sgf 902 lb_max => basis_set_b%lmax 903 lb_min => basis_set_b%lmin 904 npgfb => basis_set_b%npgf 905 nsetb = basis_set_b%nset 906 nsgfb => basis_set_b%nsgf_set 907 rpgfb => basis_set_b%pgf_radius 908 set_radius_b => basis_set_b%set_radius 909 sphi_b => basis_set_b%sphi 910 zetb => basis_set_b%zet 911 912 IF (inode == 1) THEN 913 last_jatom = 0 914 alpha = 1.0_dp 915 END IF 916 ldsa = SIZE(sphi_a, 1) 917 ldsb = SIZE(sphi_b, 1) 918 919 IF (my_minimum_image) THEN 920 IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE 921 END IF 922 923 IF (jatom /= last_jatom) THEN 924 new_atom_b = .TRUE. 925 last_jatom = jatom 926 ELSE 927 new_atom_b = .FALSE. 928 END IF 929 930 IF (new_atom_b) THEN 931 IF (iatom <= jatom) THEN 932 irow = iatom 933 icol = jatom 934 alpha = 1.0_dp 935 ELSE 936 irow = jatom 937 icol = iatom 938 IF (dbcsr_get_matrix_type(op(1)%matrix) == dbcsr_type_antisymmetric) THEN 939 !IF(op(1)%matrix%symmetry=="antisymmetric") THEN 940 alpha = -1.0_dp 941 END IF 942 END IF 943 944 DO i = 1, 3 945 NULLIFY (op_dip(i)%block) 946 CALL dbcsr_get_block_p(matrix=op(i)%matrix, & 947 row=irow, col=icol, block=op_dip(i)%block, found=found) 948 CPASSERT(ASSOCIATED(op_dip(i)%block)) 949 END DO 950 END IF ! new_atom_b 951 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 952 dab = SQRT(rab2) 953 954 DO iset = 1, nseta 955 956 ncoa = npgfa(iset)*ncoset(la_max(iset)) 957 sgfa = first_sgfa(1, iset) 958 959 DO jset = 1, nsetb 960 961 ncob = npgfb(jset)*ncoset(lb_max(jset)) 962 sgfb = first_sgfb(1, jset) 963 964 IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN 965 966! *** Calculate the primitive overlap integrals *** 967 CALL diffop(la_max(iset), npgfa(iset), zeta(:, iset), & 968 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), & 969 zetb(:, jset), rpgfb(:, jset), lb_min(jset), rab, difab) 970 971! *** Contraction *** 972 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 973 alpha, difab(1, 1, 1), ldab, sphi_b(1, sgfb), ldsb, & 974 0.0_dp, work(1, 1), ldwork) 975 IF (iatom <= jatom) THEN 976 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 977 1.0_dp, sphi_a(1, sgfa), ldsa, & 978 work(1, 1), ldwork, & 979 1.0_dp, op_dip(1)%block(sgfa, sgfb), & 980 SIZE(op_dip(1)%block, 1)) 981 982 ELSE 983 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, & 984 1.0_dp, work(1, 1), ldwork, & 985 sphi_a(1, sgfa), ldsa, & 986 1.0_dp, op_dip(1)%block(sgfb, sgfa), & 987 SIZE(op_dip(1)%block, 1)) 988 989 END IF 990 991! *** Contraction *** 992 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 993 alpha, difab(1, 1, 2), ldab, sphi_b(1, sgfb), ldsb, & 994 0.0_dp, work(1, 1), ldwork) 995 IF (iatom <= jatom) THEN 996 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 997 1.0_dp, sphi_a(1, sgfa), ldsa, & 998 work(1, 1), ldwork, & 999 1.0_dp, op_dip(2)%block(sgfa, sgfb), & 1000 SIZE(op_dip(2)%block, 1)) 1001 ELSE 1002 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, & 1003 1.0_dp, work(1, 1), ldwork, & 1004 sphi_a(1, sgfa), ldsa, & 1005 1.0_dp, op_dip(2)%block(sgfb, sgfa), & 1006 SIZE(op_dip(2)%block, 1)) 1007 END IF 1008 1009! *** Contraction *** 1010 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 1011 alpha, difab(1, 1, 3), ldab, sphi_b(1, sgfb), ldsb, & 1012 0.0_dp, work(1, 1), ldwork) 1013 IF (iatom <= jatom) THEN 1014 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 1015 1.0_dp, sphi_a(1, sgfa), ldsa, & 1016 work(1, 1), ldwork, & 1017 1.0_dp, op_dip(3)%block(sgfa, sgfb), & 1018 SIZE(op_dip(3)%block, 1)) 1019 ELSE 1020 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, & 1021 1.0_dp, work(1, 1), ldwork, & 1022 sphi_a(1, sgfa), ldsa, & 1023 1.0_dp, op_dip(3)%block(sgfb, sgfa), & 1024 SIZE(op_dip(3)%block, 1)) 1025 END IF 1026 END IF ! >= dab 1027 1028 END DO ! jset 1029 1030 END DO ! iset 1031 1032 END DO 1033 CALL neighbor_list_iterator_release(nl_iterator) 1034 1035 DO i = 1, 3 1036 NULLIFY (op_dip(i)%block) 1037 END DO 1038 DEALLOCATE (op_dip) 1039 1040 DEALLOCATE (difab, work, basis_set_list) 1041 1042 CALL timestop(handle) 1043 1044 END SUBROUTINE p_xyz_ao 1045 1046! ************************************************************************************************** 1047!> \brief Calculation of the components of the dipole operator in the length form 1048!> by taking the relative position operator r-Rc, with respect a reference point Rc 1049!> Probably it does not work for PBC, or maybe yes if the wfn are 1050!> sufficiently localized 1051!> The elements of the sparse matrices are the integrals in the 1052!> basis functions 1053!> \param op matrix representation of the p operator 1054!> calculated in terms of the contracted basis functions 1055!> \param qs_env environment for the lists and the basis sets 1056!> \param rc reference vector position 1057!> \param order maximum order of the momentum, for the dipole order = 1, order = -2 for quad only 1058!> \param minimum_image take into account only the first neighbors in the lists 1059!> \param soft ... 1060!> \par History 1061!> 03.2006 created [MI] 1062!> 06.2019 added quarupole only option (A.Bussy) 1063!> \author MI 1064! ************************************************************************************************** 1065 1066 SUBROUTINE rRc_xyz_ao(op, qs_env, rc, order, minimum_image, soft) 1067 1068 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op 1069 TYPE(qs_environment_type), POINTER :: qs_env 1070 REAL(dp) :: Rc(3) 1071 INTEGER, INTENT(IN) :: order 1072 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft 1073 1074 CHARACTER(len=*), PARAMETER :: routineN = 'rRc_xyz_ao', routineP = moduleN//':'//routineN 1075 1076 INTEGER :: handle, iatom, icol, ikind, imom, inode, irow, iset, jatom, jkind, jset, & 1077 last_jatom, ldab, ldsa, ldsb, ldwork, M_dim, maxl, ncoa, ncob, nkind, nseta, nsetb, sgfa, & 1078 sgfb, smom 1079 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, npgfa, npgfb, & 1080 nsgfa, nsgfb 1081 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 1082 LOGICAL :: found, my_minimum_image, my_soft, & 1083 new_atom_b 1084 REAL(KIND=dp) :: dab, Lxo2, Lyo2, Lzo2, rab2 1085 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc 1086 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 1087 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, & 1088 zeta, zetb 1089 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab 1090 TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip 1091 TYPE(cell_type), POINTER :: cell 1092 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 1093 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 1094 TYPE(neighbor_list_iterator_p_type), & 1095 DIMENSION(:), POINTER :: nl_iterator 1096 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1097 POINTER :: sab_orb 1098 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1099 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1100 TYPE(qs_kind_type), POINTER :: qs_kind 1101 1102 CALL timeset(routineN, handle) 1103 1104 NULLIFY (qs_kind, qs_kind_set) 1105 NULLIFY (cell, particle_set) 1106 NULLIFY (sab_orb) 1107 NULLIFY (mab, op_dip, work) 1108 NULLIFY (la_max, la_min, lb_max, npgfa, npgfb, nsgfa, nsgfb) 1109 NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb) 1110 1111 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 1112 cell=cell, particle_set=particle_set, sab_orb=sab_orb) 1113 1114 nkind = SIZE(qs_kind_set) 1115 1116 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 1117 maxco=ldwork, maxlgto=maxl) 1118 1119 my_minimum_image = .FALSE. 1120 IF (PRESENT(minimum_image)) THEN 1121 my_minimum_image = minimum_image 1122 Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp 1123 Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp 1124 Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp 1125 END IF 1126 my_soft = .FALSE. 1127 IF (PRESENT(soft)) THEN 1128 my_soft = soft 1129 END IF 1130 1131 ldab = ldwork 1132 1133 smom = 1 1134 IF (order == -2) smom = 4 1135 M_dim = ncoset(ABS(order)) - 1 1136 CPASSERT(M_dim <= SIZE(op, 1)) 1137 1138 ALLOCATE (mab(ldab, ldab, 1:M_dim)) 1139 mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp 1140 ALLOCATE (work(ldwork, ldwork)) 1141 work(1:ldwork, 1:ldwork) = 0.0_dp 1142 ALLOCATE (op_dip(smom:M_dim)) 1143 1144 DO imom = smom, M_dim 1145 NULLIFY (op_dip(imom)%block) 1146 END DO 1147 1148 ALLOCATE (basis_set_list(nkind)) 1149 DO ikind = 1, nkind 1150 qs_kind => qs_kind_set(ikind) 1151 CALL get_qs_kind(qs_kind=qs_kind, softb=my_soft, basis_set=basis_set_a) 1152 IF (ASSOCIATED(basis_set_a)) THEN 1153 basis_set_list(ikind)%gto_basis_set => basis_set_a 1154 ELSE 1155 NULLIFY (basis_set_list(ikind)%gto_basis_set) 1156 END IF 1157 END DO 1158 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 1159 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 1160 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 1161 iatom=iatom, jatom=jatom, r=rab) 1162 basis_set_a => basis_set_list(ikind)%gto_basis_set 1163 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 1164 basis_set_b => basis_set_list(jkind)%gto_basis_set 1165 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 1166 ra = pbc(particle_set(iatom)%r, cell) 1167 ! basis ikind 1168 first_sgfa => basis_set_a%first_sgf 1169 la_max => basis_set_a%lmax 1170 la_min => basis_set_a%lmin 1171 npgfa => basis_set_a%npgf 1172 nseta = basis_set_a%nset 1173 nsgfa => basis_set_a%nsgf_set 1174 rpgfa => basis_set_a%pgf_radius 1175 set_radius_a => basis_set_a%set_radius 1176 sphi_a => basis_set_a%sphi 1177 zeta => basis_set_a%zet 1178 ! basis jkind 1179 first_sgfb => basis_set_b%first_sgf 1180 lb_max => basis_set_b%lmax 1181 npgfb => basis_set_b%npgf 1182 nsetb = basis_set_b%nset 1183 nsgfb => basis_set_b%nsgf_set 1184 rpgfb => basis_set_b%pgf_radius 1185 set_radius_b => basis_set_b%set_radius 1186 sphi_b => basis_set_b%sphi 1187 zetb => basis_set_b%zet 1188 1189 ldsa = SIZE(sphi_a, 1) 1190 ldsb = SIZE(sphi_b, 1) 1191 IF (inode == 1) last_jatom = 0 1192 1193 IF (my_minimum_image) THEN 1194 IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE 1195 END IF 1196 1197 rb = rab + ra 1198 1199 IF (jatom /= last_jatom) THEN 1200 new_atom_b = .TRUE. 1201 last_jatom = jatom 1202 ELSE 1203 new_atom_b = .FALSE. 1204 END IF 1205 1206 IF (new_atom_b) THEN 1207 IF (iatom <= jatom) THEN 1208 irow = iatom 1209 icol = jatom 1210 ELSE 1211 irow = jatom 1212 icol = iatom 1213 END IF 1214 1215 DO imom = smom, M_dim 1216 NULLIFY (op_dip(imom)%block) 1217 CALL dbcsr_get_block_p(matrix=op(imom)%matrix, & 1218 row=irow, col=icol, block=op_dip(imom)%block, found=found) 1219 CPASSERT(ASSOCIATED(op_dip(imom)%block)) 1220 END DO ! imom 1221 END IF ! new_atom_b 1222 1223 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 1224 dab = SQRT(rab2) 1225 1226 DO iset = 1, nseta 1227 1228 ncoa = npgfa(iset)*ncoset(la_max(iset)) 1229 sgfa = first_sgfa(1, iset) 1230 1231 DO jset = 1, nsetb 1232 1233 ncob = npgfb(jset)*ncoset(lb_max(jset)) 1234 sgfb = first_sgfb(1, jset) 1235 1236 IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN 1237 1238 rac = pbc(rc, ra, cell) 1239 rbc = pbc(rc, rb, cell) 1240 1241! *** Calculate the primitive overlap integrals *** 1242 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), & 1243 rpgfa(:, iset), la_min(iset), & 1244 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), & 1245 ABS(order), rac, rbc, mab) 1246 1247 DO imom = smom, M_dim 1248! *** Contraction *** 1249 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 1250 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, & 1251 0.0_dp, work(1, 1), ldwork) 1252 IF (iatom <= jatom) THEN 1253 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 1254 1.0_dp, sphi_a(1, sgfa), ldsa, & 1255 work(1, 1), ldwork, & 1256 1.0_dp, op_dip(imom)%block(sgfa, sgfb), & 1257 SIZE(op_dip(imom)%block, 1)) 1258 ELSE 1259 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, & 1260 1.0_dp, work(1, 1), ldwork, & 1261 sphi_a(1, sgfa), ldsa, & 1262 1.0_dp, op_dip(imom)%block(sgfb, sgfa), & 1263 SIZE(op_dip(imom)%block, 1)) 1264 END IF 1265 1266 END DO ! imom 1267 END IF ! >= dab 1268 1269 END DO ! jset 1270 1271 END DO ! iset 1272 1273 END DO 1274 CALL neighbor_list_iterator_release(nl_iterator) 1275 1276 DO imom = smom, M_dim 1277 NULLIFY (op_dip(imom)%block) 1278 END DO 1279 DEALLOCATE (op_dip) 1280 1281 DEALLOCATE (mab, work, basis_set_list) 1282 1283 CALL timestop(handle) 1284 1285 END SUBROUTINE rRc_xyz_ao 1286 1287! ************************************************************************************************** 1288!> \brief Calculation of the multipole operators integrals 1289!> and of its derivatives of the type 1290!> [\mu | op | d(\nu)/dR(\nu)]-[d(\mu)/dR(\mu)| op | \nu] 1291!> by taking the relative position operator r-Rc, with respect a reference point Rc 1292!> The derivative are with respect to the primitive position, 1293!> The multipole operator is symmetric and if it does not depend on R(\mu) or R(\nu) 1294!> therefore [\mu | op | d(\nu)/dR(\nu)] = -[d(\mu)/dR(\mu)| op | \nu] 1295!> [\mu|op|d(\nu)/dR]-[d(\mu)/dR|op|\nu]=2[\mu|op|d(\nu)/dR] 1296!> When it is not the case a correction term is needed 1297!> 1298!> The momentum operator [\mu|M|\nu] is symmetric, the number of components is 1299!> determined by the order: 3 for order 1 (x,y,x), 9 for order 2(xx,xy,xz,yy,yz,zz) 1300!> The derivative of the type [\mu | op | d(\nu)/dR_i(\nu)], where 1301!> i indicates the cartesian direction, is antisymmetric only when 1302!> the no component M =(r_i) or (r_i r_j) is in the same cartesian 1303!> direction of the derivative, indeed 1304!> d([\mu|M|\nu])/dr_i = [d(\mu)/dr_i|M|\nu] + [\mu|M|d(\nu)/dr_i] + [\mu |d(M)/dr_i|\nu] 1305!> d([\mu|M|\nu])/dr_i = -[d(\mu)/dR_i(\mu)|M|\nu] -[\mu|M|d(\nu)/dR_i(\nu)] + [\mu |d(M)/dr_i|\nu] 1306!> Therefore we cannot use an antisymmetric matrix 1307!> 1308!> The same holds for the derivative with respect to the electronic position r 1309!> taking into account that [\mu|op|d(\nu)/dR] = -[\mu|op|d(\nu)/dr] 1310!> \param op matrix representation of the p operator 1311!> calculated in terms of the contracted basis functions 1312!> \param op_der ... 1313!> \param qs_env environment for the lists and the basis sets 1314!> \param rc reference vector position 1315!> \param order maximum order of the momentum, for the dipole order = 1 1316!> \param minimum_image take into account only the first neighbors in the lists 1317!> \param soft ... 1318!> \par History 1319!> 03.2006 created [MI] 1320!> \author MI 1321!> \note 1322!> Probably it does not work for PBC, or maybe yes if the wfn are 1323!> sufficiently localized 1324!> The elements of the sparse matrices are the integrals in the 1325!> basis functions 1326! ************************************************************************************************** 1327 SUBROUTINE rRc_xyz_der_ao(op, op_der, qs_env, rc, order, minimum_image, soft) 1328 1329 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op 1330 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_der 1331 TYPE(qs_environment_type), POINTER :: qs_env 1332 REAL(dp) :: Rc(3) 1333 INTEGER, INTENT(IN) :: order 1334 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft 1335 1336 CHARACTER(len=*), PARAMETER :: routineN = 'rRc_xyz_der_ao', routineP = moduleN//':'//routineN 1337 1338 INTEGER :: handle, i, iatom, icol, idir, ikind, imom, inode, ipgf, irow, iset, j, jatom, & 1339 jkind, jpgf, jset, last_jatom, lda_min, ldab, ldb_min, ldsa, ldsb, ldwork, M_dim, maxl, & 1340 na, nb, ncoa, ncob, nda, ndb, nkind, nseta, nsetb, sgfa, sgfb 1341 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 1342 npgfb, nsgfa, nsgfb 1343 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 1344 LOGICAL :: my_minimum_image, my_soft, new_atom_b, & 1345 op_der_found, op_found 1346 REAL(KIND=dp) :: alpha, alpha_der, dab, Lxo2, Lyo2, Lzo2, & 1347 rab2 1348 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc 1349 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 1350 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, & 1351 zeta, zetb 1352 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab, mab_tmp 1353 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: difmab 1354 TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip 1355 TYPE(block_p_type), DIMENSION(:, :), POINTER :: op_dip_der 1356 TYPE(cell_type), POINTER :: cell 1357 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 1358 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 1359 TYPE(neighbor_list_iterator_p_type), & 1360 DIMENSION(:), POINTER :: nl_iterator 1361 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1362 POINTER :: sab_all 1363 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1364 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1365 TYPE(qs_kind_type), POINTER :: qs_kind 1366 1367 CALL timeset(routineN, handle) 1368 1369 CPASSERT(ASSOCIATED(op)) 1370 CPASSERT(ASSOCIATED(op_der)) 1371 !IF(.NOT.op_sm_der(1,1)%matrix%symmetry=="none") THEN 1372 IF (.NOT. dbcsr_get_matrix_type(op_der(1, 1)%matrix) .EQ. dbcsr_type_no_symmetry) THEN 1373 CPABORT("") 1374 END IF 1375 1376 NULLIFY (qs_kind, qs_kind_set) 1377 NULLIFY (cell, particle_set) 1378 NULLIFY (sab_all) 1379 NULLIFY (difmab, mab, mab_tmp) 1380 NULLIFY (op_dip, op_dip_der, work) 1381 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb) 1382 NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb) 1383 1384 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 1385 cell=cell, particle_set=particle_set, & 1386 sab_all=sab_all) 1387 1388 nkind = SIZE(qs_kind_set) 1389 1390 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 1391 maxco=ldwork, maxlgto=maxl) 1392 1393 my_minimum_image = .FALSE. 1394 IF (PRESENT(minimum_image)) THEN 1395 my_minimum_image = minimum_image 1396 Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp 1397 Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp 1398 Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp 1399 END IF 1400 my_soft = .FALSE. 1401 IF (PRESENT(soft)) THEN 1402 my_soft = soft 1403 END IF 1404 1405 ldab = ldwork 1406 1407 M_dim = ncoset(order) - 1 1408 CPASSERT(M_dim <= SIZE(op, 1)) 1409 1410 ALLOCATE (mab(ldab, ldab, M_dim)) 1411 mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp 1412 ALLOCATE (difmab(ldab, ldab, M_dim, 3)) 1413 difmab(1:ldab, 1:ldab, 1:M_dim, 1:3) = 0.0_dp 1414 1415 ALLOCATE (work(ldwork, ldwork)) 1416 work(1:ldwork, 1:ldwork) = 0.0_dp 1417 ALLOCATE (op_dip(M_dim)) 1418 ALLOCATE (op_dip_der(M_dim, 3)) 1419 1420 DO imom = 1, M_dim 1421 NULLIFY (op_dip(imom)%block) 1422 DO i = 1, 3 1423 NULLIFY (op_dip_der(imom, i)%block) 1424 END DO 1425 END DO 1426 1427 ALLOCATE (basis_set_list(nkind)) 1428 DO ikind = 1, nkind 1429 qs_kind => qs_kind_set(ikind) 1430 CALL get_qs_kind(qs_kind=qs_kind, softb=my_soft, basis_set=basis_set_a) 1431 IF (ASSOCIATED(basis_set_a)) THEN 1432 basis_set_list(ikind)%gto_basis_set => basis_set_a 1433 ELSE 1434 NULLIFY (basis_set_list(ikind)%gto_basis_set) 1435 END IF 1436 END DO 1437 CALL neighbor_list_iterator_create(nl_iterator, sab_all) 1438 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 1439 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 1440 iatom=iatom, jatom=jatom, r=rab) 1441 basis_set_a => basis_set_list(ikind)%gto_basis_set 1442 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 1443 basis_set_b => basis_set_list(jkind)%gto_basis_set 1444 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 1445 ra = pbc(particle_set(iatom)%r, cell) 1446 ! basis ikind 1447 first_sgfa => basis_set_a%first_sgf 1448 la_max => basis_set_a%lmax 1449 la_min => basis_set_a%lmin 1450 npgfa => basis_set_a%npgf 1451 nseta = basis_set_a%nset 1452 nsgfa => basis_set_a%nsgf_set 1453 rpgfa => basis_set_a%pgf_radius 1454 set_radius_a => basis_set_a%set_radius 1455 sphi_a => basis_set_a%sphi 1456 zeta => basis_set_a%zet 1457 ! basis jkind 1458 first_sgfb => basis_set_b%first_sgf 1459 lb_max => basis_set_b%lmax 1460 lb_min => basis_set_b%lmin 1461 npgfb => basis_set_b%npgf 1462 nsetb = basis_set_b%nset 1463 nsgfb => basis_set_b%nsgf_set 1464 rpgfb => basis_set_b%pgf_radius 1465 set_radius_b => basis_set_b%set_radius 1466 sphi_b => basis_set_b%sphi 1467 zetb => basis_set_b%zet 1468 1469 ldsa = SIZE(sphi_a, 1) 1470 IF (ldsa .EQ. 0) CYCLE 1471 ldsb = SIZE(sphi_b, 1) 1472 IF (ldsb .EQ. 0) CYCLE 1473 IF (inode == 1) last_jatom = 0 1474 1475 IF (my_minimum_image) THEN 1476 IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE 1477 END IF 1478 1479 rb = rab + ra 1480 1481 IF (jatom /= last_jatom) THEN 1482 new_atom_b = .TRUE. 1483 last_jatom = jatom 1484 ELSE 1485 new_atom_b = .FALSE. 1486 END IF 1487 1488 IF (new_atom_b) THEN 1489 irow = iatom 1490 icol = jatom 1491 alpha_der = 2.0_dp 1492 1493 DO imom = 1, M_dim 1494 NULLIFY (op_dip(imom)%block) 1495 CALL dbcsr_get_block_p(matrix=op(imom)%matrix, & 1496 row=irow, col=icol, & 1497 block=op_dip(imom)%block, & 1498 found=op_found) 1499 CPASSERT(op_found .AND. ASSOCIATED(op_dip(imom)%block)) 1500 DO idir = 1, 3 1501 NULLIFY (op_dip_der(imom, idir)%block) 1502 CALL dbcsr_get_block_p(matrix=op_der(imom, idir)%matrix, & 1503 row=irow, col=icol, & 1504 block=op_dip_der(imom, idir)%block, & 1505 found=op_der_found) 1506 CPASSERT(op_der_found .AND. ASSOCIATED(op_dip_der(imom, idir)%block)) 1507 END DO ! idir 1508 END DO ! imom 1509 END IF ! new_atom_b 1510 1511 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 1512 dab = SQRT(rab2) 1513 1514 DO iset = 1, nseta 1515 1516 ncoa = npgfa(iset)*ncoset(la_max(iset)) 1517 sgfa = first_sgfa(1, iset) 1518 1519 DO jset = 1, nsetb 1520 1521 ncob = npgfb(jset)*ncoset(lb_max(jset)) 1522 sgfb = first_sgfb(1, jset) 1523 1524 IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN 1525 1526 rac = pbc(rc, ra, cell) 1527 rbc = rac + rab 1528! rac = pbc(rc,ra,cell) 1529! rbc = pbc(rc,rb,cell) 1530 1531 ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), & 1532 npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(order) - 1)) 1533 1534 lda_min = MAX(0, la_min(iset) - 1) 1535 ldb_min = MAX(0, lb_min(jset) - 1) 1536! *** Calculate the primitive overlap integrals *** 1537 CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), & 1538 rpgfa(:, iset), lda_min, & 1539 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), rpgfb(:, jset), & 1540 order, rac, rbc, mab_tmp) 1541 1542! *** Calculate the derivatives 1543 CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), & 1544 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), & 1545 zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, & 1546 difmab, mab_ext=mab_tmp) 1547 1548! Contract and copy in the sparse matrix 1549 mab = 0.0_dp 1550 DO imom = 1, M_dim 1551 na = 0 1552 nda = 0 1553 DO ipgf = 1, npgfa(iset) 1554 nb = 0 1555 ndb = 0 1556 DO jpgf = 1, npgfb(jset) 1557 DO j = 1, ncoset(lb_max(jset)) 1558 DO i = 1, ncoset(la_max(iset)) 1559 mab(i + na, j + nb, imom) = mab_tmp(i + nda, j + ndb, imom) 1560 END DO ! i 1561 END DO ! j 1562 nb = nb + ncoset(lb_max(jset)) 1563 ndb = ndb + ncoset(lb_max(jset) + 1) 1564 END DO ! jpgf 1565 na = na + ncoset(la_max(iset)) 1566 nda = nda + ncoset(la_max(iset) + 1) 1567 END DO ! ipgf 1568 1569! *** Contraction *** 1570 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 1571 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, & 1572 0.0_dp, work(1, 1), ldwork) 1573 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 1574 1.0_dp, sphi_a(1, sgfa), ldsa, & 1575 work(1, 1), ldwork, & 1576 1.0_dp, op_dip(imom)%block(sgfa, sgfb), & 1577 SIZE(op_dip(imom)%block, 1)) 1578 1579 alpha = -1.0_dp !-alpha_der 1580 DO idir = 1, 3 1581 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 1582 alpha, difmab(1, 1, imom, idir), ldab, sphi_b(1, sgfb), ldsb, & 1583 0.0_dp, work(1, 1), ldwork) 1584 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 1585 1.0_dp, sphi_a(1, sgfa), ldsa, & 1586 work(1, 1), ldwork, & 1587 1.0_dp, op_dip_der(imom, idir)%block(sgfa, sgfb), & 1588 SIZE(op_dip_der(imom, idir)%block, 1)) 1589 1590 END DO ! idir 1591 1592 END DO ! imom 1593 1594 DEALLOCATE (mab_tmp) 1595 END IF ! >= dab 1596 1597 END DO ! jset 1598 1599 END DO ! iset 1600 1601 END DO 1602 CALL neighbor_list_iterator_release(nl_iterator) 1603 1604 DO i = 1, 3 1605 NULLIFY (op_dip(i)%block) 1606 END DO 1607 DEALLOCATE (op_dip, op_dip_der) 1608 1609 DEALLOCATE (mab, difmab, work, basis_set_list) 1610 1611 CALL timestop(handle) 1612 1613 END SUBROUTINE rRc_xyz_der_ao 1614 1615END MODULE qs_operators_ao 1616 1617