1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate the atomic operator matrices 8!> \author jgh 9!> \date 03.03.2008 10!> \version 1.0 11!> 12! ************************************************************************************************** 13MODULE atom_operators 14 USE ai_onecenter, ONLY: & 15 sg_coulomb, sg_erf, sg_erfc, sg_exchange, sg_gpot, sg_kinetic, sg_kinnuc, sg_nuclear, & 16 sg_overlap, sg_proj_ol, sto_kinetic, sto_nuclear, sto_overlap 17 USE atom_types, ONLY: & 18 atom_basis_gridrep, atom_basis_type, atom_compare_grids, atom_integrals, & 19 atom_potential_type, atom_state, cgto_basis, ecp_pseudo, gth_pseudo, gto_basis, lmat, & 20 no_pseudo, num_basis, release_atom_basis, sgp_pseudo, sto_basis, upf_pseudo 21 USE atom_utils, ONLY: & 22 atom_solve, contract2, contract2add, contract4, coulomb_potential_numeric, integrate_grid, & 23 numpot_matrix, slater_density, wigner_slater_functional 24 USE dkh_main, ONLY: dkh_atom_transformation 25 USE input_constants, ONLY: & 26 barrier_conf, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, & 27 do_sczoramp_atom, do_zoramp_atom, poly_conf 28 USE kinds, ONLY: dp 29 USE lapack, ONLY: lapack_sgesv 30 USE mathconstants, ONLY: gamma1,& 31 sqrt2 32 USE mathlib, ONLY: jacobi 33 USE periodic_table, ONLY: ptable 34 USE physcon, ONLY: c_light_au 35 USE qs_grid_atom, ONLY: grid_atom_type 36#include "./base/base_uses.f90" 37 38 IMPLICIT NONE 39 40 PRIVATE 41 42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_operators' 43 44 PUBLIC :: atom_int_setup, atom_ppint_setup, atom_int_release, atom_ppint_release 45 PUBLIC :: atom_relint_setup, atom_relint_release, atom_basis_projection_overlap 46 PUBLIC :: calculate_model_potential 47 48CONTAINS 49 50! ************************************************************************************************** 51!> \brief Set up atomic integrals. 52!> \param integrals atomic integrals 53!> \param basis atomic basis set 54!> \param potential pseudo-potential 55!> \param eri_coulomb setup one-centre Coulomb Electron Repulsion Integrals 56!> \param eri_exchange setup one-centre exchange Electron Repulsion Integrals 57!> \param all_nu compute integrals for all even integer parameters [0 .. 2*l] 58!> REDUNDANT, AS THIS SUBROUTINE IS NEVER INVOKED WITH all_nu = .TRUE. 59! ************************************************************************************************** 60 SUBROUTINE atom_int_setup(integrals, basis, potential, & 61 eri_coulomb, eri_exchange, all_nu) 62 TYPE(atom_integrals), INTENT(INOUT) :: integrals 63 TYPE(atom_basis_type), INTENT(INOUT) :: basis 64 TYPE(atom_potential_type), INTENT(IN), OPTIONAL :: potential 65 LOGICAL, INTENT(IN), OPTIONAL :: eri_coulomb, eri_exchange, all_nu 66 67 CHARACTER(len=*), PARAMETER :: routineN = 'atom_int_setup', routineP = moduleN//':'//routineN 68 69 INTEGER :: handle, i, ii, info, ipiv(1000), l, l1, & 70 l2, ll, lwork, m, m1, m2, mm1, mm2, n, & 71 n1, n2, nn1, nn2, nu, nx 72 REAL(KIND=dp) :: om, rc, ron, sc, x 73 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, w, work 74 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat, vmat 75 REAL(KIND=dp), DIMENSION(:, :), POINTER :: eri 76 77 CALL timeset(routineN, handle) 78 79 IF (integrals%status == 0) THEN 80 n = MAXVAL(basis%nbas) 81 integrals%n = basis%nbas 82 83 IF (PRESENT(eri_coulomb)) THEN 84 integrals%eri_coulomb = eri_coulomb 85 ELSE 86 integrals%eri_coulomb = .FALSE. 87 END IF 88 IF (PRESENT(eri_exchange)) THEN 89 integrals%eri_exchange = eri_exchange 90 ELSE 91 integrals%eri_exchange = .FALSE. 92 END IF 93 IF (PRESENT(all_nu)) THEN 94 integrals%all_nu = all_nu 95 ELSE 96 integrals%all_nu = .FALSE. 97 END IF 98 99 NULLIFY (integrals%ovlp, integrals%kin, integrals%core, integrals%conf) 100 DO ll = 1, SIZE(integrals%ceri) 101 NULLIFY (integrals%ceri(ll)%int, integrals%eeri(ll)%int) 102 END DO 103 104 ALLOCATE (integrals%ovlp(n, n, 0:lmat)) 105 integrals%ovlp = 0._dp 106 107 ALLOCATE (integrals%kin(n, n, 0:lmat)) 108 integrals%kin = 0._dp 109 110 integrals%status = 1 111 112 IF (PRESENT(potential)) THEN 113 IF (potential%confinement) THEN 114 ALLOCATE (integrals%conf(n, n, 0:lmat)) 115 integrals%conf = 0._dp 116 m = basis%grid%nr 117 ALLOCATE (cpot(1:m)) 118 IF (potential%conf_type == poly_conf) THEN 119 rc = potential%rcon 120 sc = potential%scon 121 cpot(1:m) = (basis%grid%rad(1:m)/rc)**sc 122 ELSEIF (potential%conf_type == barrier_conf) THEN 123 om = potential%rcon 124 ron = potential%scon 125 rc = ron + om 126 DO i = 1, m 127 IF (basis%grid%rad(i) < ron) THEN 128 cpot(i) = 0.0_dp 129 ELSEIF (basis%grid%rad(i) < rc) THEN 130 x = (basis%grid%rad(i) - ron)/om 131 x = 1._dp - x 132 cpot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp 133 x = (rc - basis%grid%rad(i))**2/om/(basis%grid%rad(i) - ron) 134 cpot(i) = cpot(i)*x 135 ELSE 136 cpot(i) = 1.0_dp 137 END IF 138 END DO 139 ELSE 140 CPABORT("") 141 END IF 142 CALL numpot_matrix(integrals%conf, cpot, basis, 0) 143 DEALLOCATE (cpot) 144 END IF 145 END IF 146 147 SELECT CASE (basis%basis_type) 148 CASE DEFAULT 149 CPABORT("") 150 CASE (GTO_BASIS) 151 DO l = 0, lmat 152 n = integrals%n(l) 153 CALL sg_overlap(integrals%ovlp(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l)) 154 CALL sg_kinetic(integrals%kin(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l)) 155 END DO 156 IF (integrals%eri_coulomb) THEN 157 ll = 0 158 DO l1 = 0, lmat 159 n1 = integrals%n(l1) 160 nn1 = (n1*(n1 + 1))/2 161 DO l2 = 0, l1 162 n2 = integrals%n(l2) 163 nn2 = (n2*(n2 + 1))/2 164 IF (integrals%all_nu) THEN 165 nx = MIN(2*l1, 2*l2) 166 ELSE 167 nx = 0 168 END IF 169 DO nu = 0, nx, 2 170 ll = ll + 1 171 CPASSERT(ll <= SIZE(integrals%ceri)) 172 ALLOCATE (integrals%ceri(ll)%int(nn1, nn2)) 173 integrals%ceri(ll)%int = 0._dp 174 eri => integrals%ceri(ll)%int 175 CALL sg_coulomb(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2) 176 END DO 177 END DO 178 END DO 179 END IF 180 IF (integrals%eri_exchange) THEN 181 ll = 0 182 DO l1 = 0, lmat 183 n1 = integrals%n(l1) 184 nn1 = (n1*(n1 + 1))/2 185 DO l2 = 0, l1 186 n2 = integrals%n(l2) 187 nn2 = (n2*(n2 + 1))/2 188 DO nu = ABS(l1 - l2), l1 + l2, 2 189 ll = ll + 1 190 CPASSERT(ll <= SIZE(integrals%eeri)) 191 ALLOCATE (integrals%eeri(ll)%int(nn1, nn2)) 192 integrals%eeri(ll)%int = 0._dp 193 eri => integrals%eeri(ll)%int 194 CALL sg_exchange(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2) 195 END DO 196 END DO 197 END DO 198 END IF 199 CASE (CGTO_BASIS) 200 DO l = 0, lmat 201 n = integrals%n(l) 202 m = basis%nprim(l) 203 IF (n > 0 .AND. m > 0) THEN 204 ALLOCATE (omat(m, m)) 205 206 CALL sg_overlap(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l)) 207 CALL contract2(integrals%ovlp(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l)) 208 CALL sg_kinetic(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l)) 209 CALL contract2(integrals%kin(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l)) 210 DEALLOCATE (omat) 211 END IF 212 END DO 213 IF (integrals%eri_coulomb) THEN 214 ll = 0 215 DO l1 = 0, lmat 216 n1 = integrals%n(l1) 217 nn1 = (n1*(n1 + 1))/2 218 m1 = basis%nprim(l1) 219 mm1 = (m1*(m1 + 1))/2 220 DO l2 = 0, l1 221 n2 = integrals%n(l2) 222 nn2 = (n2*(n2 + 1))/2 223 m2 = basis%nprim(l2) 224 mm2 = (m2*(m2 + 1))/2 225 IF (integrals%all_nu) THEN 226 nx = MIN(2*l1, 2*l2) 227 ELSE 228 nx = 0 229 END IF 230 DO nu = 0, nx, 2 231 ll = ll + 1 232 CPASSERT(ll <= SIZE(integrals%ceri)) 233 ALLOCATE (integrals%ceri(ll)%int(nn1, nn2)) 234 integrals%ceri(ll)%int = 0._dp 235 ALLOCATE (omat(mm1, mm2)) 236 237 eri => integrals%ceri(ll)%int 238 CALL sg_coulomb(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2) 239 CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2)) 240 241 DEALLOCATE (omat) 242 END DO 243 END DO 244 END DO 245 END IF 246 IF (integrals%eri_exchange) THEN 247 ll = 0 248 DO l1 = 0, lmat 249 n1 = integrals%n(l1) 250 nn1 = (n1*(n1 + 1))/2 251 m1 = basis%nprim(l1) 252 mm1 = (m1*(m1 + 1))/2 253 DO l2 = 0, l1 254 n2 = integrals%n(l2) 255 nn2 = (n2*(n2 + 1))/2 256 m2 = basis%nprim(l2) 257 mm2 = (m2*(m2 + 1))/2 258 DO nu = ABS(l1 - l2), l1 + l2, 2 259 ll = ll + 1 260 CPASSERT(ll <= SIZE(integrals%eeri)) 261 ALLOCATE (integrals%eeri(ll)%int(nn1, nn2)) 262 integrals%eeri(ll)%int = 0._dp 263 ALLOCATE (omat(mm1, mm2)) 264 265 eri => integrals%eeri(ll)%int 266 CALL sg_exchange(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2) 267 CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2)) 268 269 DEALLOCATE (omat) 270 END DO 271 END DO 272 END DO 273 END IF 274 CASE (STO_BASIS) 275 DO l = 0, lmat 276 n = integrals%n(l) 277 CALL sto_overlap(integrals%ovlp(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), & 278 basis%ns(1:n, l), basis%as(1:n, l)) 279 CALL sto_kinetic(integrals%kin(1:n, 1:n, l), l, basis%ns(1:n, l), basis%as(1:n, l), & 280 basis%ns(1:n, l), basis%as(1:n, l)) 281 END DO 282 CPASSERT(.NOT. integrals%eri_coulomb) 283 CPASSERT(.NOT. integrals%eri_exchange) 284 CASE (NUM_BASIS) 285 CPABORT("") 286 END SELECT 287 288 ! setup transformation matrix to get an orthogonal basis, remove linear dependencies 289 NULLIFY (integrals%utrans, integrals%uptrans) 290 n = MAXVAL(basis%nbas) 291 ALLOCATE (integrals%utrans(n, n, 0:lmat), integrals%uptrans(n, n, 0:lmat)) 292 integrals%utrans = 0._dp 293 integrals%uptrans = 0._dp 294 integrals%nne = integrals%n 295 lwork = 10*n 296 ALLOCATE (omat(n, n), w(n), vmat(n, n), work(lwork)) 297 DO l = 0, lmat 298 n = integrals%n(l) 299 IF (n > 0) THEN 300 omat(1:n, 1:n) = integrals%ovlp(1:n, 1:n, l) 301 CALL jacobi(omat(1:n, 1:n), w(1:n), vmat(1:n, 1:n)) 302 omat(1:n, 1:n) = vmat(1:n, 1:n) 303 ii = 0 304 DO i = 1, n 305 IF (w(i) > basis%eps_eig) THEN 306 ii = ii + 1 307 integrals%utrans(1:n, ii, l) = omat(1:n, i)/SQRT(w(i)) 308 END IF 309 END DO 310 integrals%nne(l) = ii 311 IF (ii > 0) THEN 312 omat(1:ii, 1:ii) = MATMUL(TRANSPOSE(integrals%utrans(1:n, 1:ii, l)), integrals%utrans(1:n, 1:ii, l)) 313 DO i = 1, ii 314 integrals%uptrans(i, i, l) = 1._dp 315 ENDDO 316 CALL lapack_sgesv(ii, ii, omat(1:ii, 1:ii), ii, ipiv, integrals%uptrans(1:ii, 1:ii, l), ii, info) 317 CPASSERT(info == 0) 318 END IF 319 END IF 320 END DO 321 DEALLOCATE (omat, vmat, w, work) 322 323 END IF 324 325 CALL timestop(handle) 326 327 END SUBROUTINE atom_int_setup 328 329! ************************************************************************************************** 330!> \brief ... 331!> \param integrals ... 332!> \param basis ... 333!> \param potential ... 334! ************************************************************************************************** 335 SUBROUTINE atom_ppint_setup(integrals, basis, potential) 336 TYPE(atom_integrals), INTENT(INOUT) :: integrals 337 TYPE(atom_basis_type), INTENT(INOUT) :: basis 338 TYPE(atom_potential_type), INTENT(IN) :: potential 339 340 CHARACTER(len=*), PARAMETER :: routineN = 'atom_ppint_setup', & 341 routineP = moduleN//':'//routineN 342 343 INTEGER :: handle, i, ii, j, k, l, m, n 344 REAL(KIND=dp) :: al, alpha, rc 345 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, xmat 346 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat, spmat 347 REAL(KIND=dp), DIMENSION(:), POINTER :: rad 348 349 CALL timeset(routineN, handle) 350 351 IF (integrals%ppstat == 0) THEN 352 n = MAXVAL(basis%nbas) 353 integrals%n = basis%nbas 354 355 NULLIFY (integrals%core, integrals%hnl) 356 357 ALLOCATE (integrals%hnl(n, n, 0:lmat)) 358 integrals%hnl = 0._dp 359 360 ALLOCATE (integrals%core(n, n, 0:lmat)) 361 integrals%core = 0._dp 362 363 ALLOCATE (integrals%clsd(n, n, 0:lmat)) 364 integrals%clsd = 0._dp 365 366 integrals%ppstat = 1 367 368 SELECT CASE (basis%basis_type) 369 CASE DEFAULT 370 CPABORT("") 371 CASE (GTO_BASIS) 372 373 SELECT CASE (potential%ppot_type) 374 CASE (no_pseudo, ecp_pseudo) 375 DO l = 0, lmat 376 n = integrals%n(l) 377 CALL sg_nuclear(integrals%core(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l)) 378 END DO 379 CASE (gth_pseudo) 380 alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp) 381 DO l = 0, lmat 382 n = integrals%n(l) 383 ALLOCATE (omat(n, n), spmat(n, 5)) 384 385 omat = 0._dp 386 CALL sg_erf(omat(1:n, 1:n), l, alpha, basis%am(1:n, l), basis%am(1:n, l)) 387 integrals%core(1:n, 1:n, l) = -potential%gth_pot%zion*omat(1:n, 1:n) 388 DO i = 1, potential%gth_pot%ncl 389 omat = 0._dp 390 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%rc, l, basis%am(1:n, l), basis%am(1:n, l)) 391 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + & 392 potential%gth_pot%cl(i)*omat(1:n, 1:n) 393 END DO 394 IF (potential%gth_pot%lpotextended) THEN 395 DO k = 1, potential%gth_pot%nexp_lpot 396 DO i = 1, potential%gth_pot%nct_lpot(k) 397 omat = 0._dp 398 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lpot(k), l, & 399 basis%am(1:n, l), basis%am(1:n, l)) 400 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + & 401 potential%gth_pot%cval_lpot(i, k)*omat(1:n, 1:n) 402 END DO 403 END DO 404 END IF 405 IF (potential%gth_pot%lsdpot) THEN 406 DO k = 1, potential%gth_pot%nexp_lsd 407 DO i = 1, potential%gth_pot%nct_lsd(k) 408 omat = 0._dp 409 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lsd(k), l, & 410 basis%am(1:n, l), basis%am(1:n, l)) 411 integrals%clsd(1:n, 1:n, l) = integrals%clsd(1:n, 1:n, l) + & 412 potential%gth_pot%cval_lsd(i, k)*omat(1:n, 1:n) 413 END DO 414 END DO 415 END IF 416 417 spmat = 0._dp 418 m = potential%gth_pot%nl(l) 419 DO i = 1, m 420 CALL sg_proj_ol(spmat(1:n, i), l, basis%am(1:n, l), i - 1, potential%gth_pot%rcnl(l)) 421 END DO 422 integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:m), & 423 MATMUL(potential%gth_pot%hnl(1:m, 1:m, l), TRANSPOSE(spmat(1:n, 1:m)))) 424 425 DEALLOCATE (omat, spmat) 426 END DO 427 CASE (upf_pseudo) 428 CALL upfint_setup(integrals, basis, potential) 429 CASE (sgp_pseudo) 430 CALL sgpint_setup(integrals, basis, potential) 431 CASE DEFAULT 432 CPABORT("") 433 END SELECT 434 435 CASE (CGTO_BASIS) 436 437 SELECT CASE (potential%ppot_type) 438 CASE (no_pseudo, ecp_pseudo) 439 DO l = 0, lmat 440 n = integrals%n(l) 441 m = basis%nprim(l) 442 ALLOCATE (omat(m, m)) 443 444 CALL sg_nuclear(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l)) 445 CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l)) 446 447 DEALLOCATE (omat) 448 END DO 449 CASE (gth_pseudo) 450 alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp) 451 DO l = 0, lmat 452 n = integrals%n(l) 453 m = basis%nprim(l) 454 IF (n > 0 .AND. m > 0) THEN 455 ALLOCATE (omat(m, m), spmat(n, 5), xmat(m)) 456 457 omat = 0._dp 458 CALL sg_erf(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l)) 459 omat(1:m, 1:m) = -potential%gth_pot%zion*omat(1:m, 1:m) 460 CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l)) 461 DO i = 1, potential%gth_pot%ncl 462 omat = 0._dp 463 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%rc, l, basis%am(1:m, l), basis%am(1:m, l)) 464 omat(1:m, 1:m) = potential%gth_pot%cl(i)*omat(1:m, 1:m) 465 CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l)) 466 END DO 467 IF (potential%gth_pot%lpotextended) THEN 468 DO k = 1, potential%gth_pot%nexp_lpot 469 DO i = 1, potential%gth_pot%nct_lpot(k) 470 omat = 0._dp 471 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lpot(k), l, & 472 basis%am(1:m, l), basis%am(1:m, l)) 473 omat(1:m, 1:m) = potential%gth_pot%cval_lpot(i, k)*omat(1:m, 1:m) 474 CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l)) 475 END DO 476 END DO 477 END IF 478 IF (potential%gth_pot%lsdpot) THEN 479 DO k = 1, potential%gth_pot%nexp_lsd 480 DO i = 1, potential%gth_pot%nct_lsd(k) 481 omat = 0._dp 482 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lsd(k), l, & 483 basis%am(1:m, l), basis%am(1:m, l)) 484 omat(1:m, 1:m) = potential%gth_pot%cval_lsd(i, k)*omat(1:m, 1:m) 485 CALL contract2add(integrals%clsd(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l)) 486 END DO 487 END DO 488 END IF 489 490 spmat = 0._dp 491 k = potential%gth_pot%nl(l) 492 DO i = 1, k 493 CALL sg_proj_ol(xmat(1:m), l, basis%am(1:m, l), i - 1, potential%gth_pot%rcnl(l)) 494 spmat(1:n, i) = MATMUL(TRANSPOSE(basis%cm(1:m, 1:n, l)), xmat(1:m)) 495 END DO 496 IF (k > 0) THEN 497 integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:k), & 498 MATMUL(potential%gth_pot%hnl(1:k, 1:k, l), & 499 TRANSPOSE(spmat(1:n, 1:k)))) 500 END IF 501 DEALLOCATE (omat, spmat, xmat) 502 END IF 503 END DO 504 CASE (upf_pseudo) 505 CALL upfint_setup(integrals, basis, potential) 506 CASE (sgp_pseudo) 507 CALL sgpint_setup(integrals, basis, potential) 508 CASE DEFAULT 509 CPABORT("") 510 END SELECT 511 512 CASE (STO_BASIS) 513 514 SELECT CASE (potential%ppot_type) 515 CASE (no_pseudo, ecp_pseudo) 516 DO l = 0, lmat 517 n = integrals%n(l) 518 CALL sto_nuclear(integrals%core(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), & 519 basis%ns(1:n, l), basis%as(1:n, l)) 520 END DO 521 CASE (gth_pseudo) 522 rad => basis%grid%rad 523 m = basis%grid%nr 524 ALLOCATE (cpot(1:m)) 525 rc = potential%gth_pot%rc 526 alpha = 1._dp/rc/SQRT(2._dp) 527 528 ! local pseudopotential, we use erf = 1/r - erfc 529 integrals%core = 0._dp 530 DO i = 1, m 531 cpot(i) = potential%gth_pot%zion*erfc(alpha*rad(i))/rad(i) 532 END DO 533 DO i = 1, potential%gth_pot%ncl 534 ii = 2*(i - 1) 535 cpot(1:m) = cpot(1:m) + potential%gth_pot%cl(i)*(rad/rc)**ii*EXP(-0.5_dp*(rad/rc)**2) 536 END DO 537 IF (potential%gth_pot%lpotextended) THEN 538 DO k = 1, potential%gth_pot%nexp_lpot 539 al = potential%gth_pot%alpha_lpot(k) 540 DO i = 1, potential%gth_pot%nct_lpot(k) 541 ii = 2*(i - 1) 542 cpot(1:m) = cpot(1:m) + potential%gth_pot%cval_lpot(i, k)*(rad/al)**ii*EXP(-0.5_dp*(rad/al)**2) 543 END DO 544 END DO 545 END IF 546 CALL numpot_matrix(integrals%core, cpot, basis, 0) 547 DO l = 0, lmat 548 n = integrals%n(l) 549 ALLOCATE (omat(n, n)) 550 omat = 0._dp 551 CALL sto_nuclear(omat(1:n, 1:n), basis%ns(1:n, l), basis%as(1:n, l), & 552 basis%ns(1:n, l), basis%as(1:n, l)) 553 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) - potential%gth_pot%zion*omat(1:n, 1:n) 554 DEALLOCATE (omat) 555 END DO 556 557 IF (potential%gth_pot%lsdpot) THEN 558 cpot = 0._dp 559 DO k = 1, potential%gth_pot%nexp_lsd 560 al = potential%gth_pot%alpha_lsd(k) 561 DO i = 1, potential%gth_pot%nct_lsd(k) 562 ii = 2*(i - 1) 563 cpot(:) = cpot + potential%gth_pot%cval_lsd(i, k)*(rad/al)**ii*EXP(-0.5_dp*(rad/al)**2) 564 END DO 565 END DO 566 CALL numpot_matrix(integrals%clsd, cpot, basis, 0) 567 END IF 568 569 DO l = 0, lmat 570 n = integrals%n(l) 571 ! non local pseudopotential 572 ALLOCATE (spmat(n, 5)) 573 spmat = 0._dp 574 k = potential%gth_pot%nl(l) 575 DO i = 1, k 576 rc = potential%gth_pot%rcnl(l) 577 cpot(:) = sqrt2/SQRT(gamma1(l + 2*i - 1))*rad**(l + 2*i - 2)*EXP(-0.5_dp*(rad/rc)**2)/rc**(l + 2*i - 0.5_dp) 578 DO j = 1, basis%nbas(l) 579 spmat(j, i) = integrate_grid(cpot, basis%bf(:, j, l), basis%grid) 580 END DO 581 END DO 582 integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:k), & 583 MATMUL(potential%gth_pot%hnl(1:k, 1:k, l), & 584 TRANSPOSE(spmat(1:n, 1:k)))) 585 DEALLOCATE (spmat) 586 END DO 587 DEALLOCATE (cpot) 588 589 CASE (upf_pseudo) 590 CALL upfint_setup(integrals, basis, potential) 591 CASE (sgp_pseudo) 592 CALL sgpint_setup(integrals, basis, potential) 593 CASE DEFAULT 594 CPABORT("") 595 END SELECT 596 597 CASE (NUM_BASIS) 598 CPABORT("") 599 END SELECT 600 601 ! add ecp_pseudo using numerical representation of basis 602 IF (potential%ppot_type == ecp_pseudo) THEN 603 ! scale 1/r potential 604 integrals%core = -potential%ecp_pot%zion*integrals%core 605 ! local potential 606 m = basis%grid%nr 607 rad => basis%grid%rad 608 ALLOCATE (cpot(1:m)) 609 cpot = 0._dp 610 DO k = 1, potential%ecp_pot%nloc 611 n = potential%ecp_pot%nrloc(k) 612 alpha = potential%ecp_pot%bloc(k) 613 cpot(:) = cpot + potential%ecp_pot%aloc(k)*rad**(n - 2)*EXP(-alpha*rad**2) 614 END DO 615 CALL numpot_matrix(integrals%core, cpot, basis, 0) 616 ! non local pseudopotential 617 DO l = 0, MIN(potential%ecp_pot%lmax, lmat) 618 cpot = 0._dp 619 DO k = 1, potential%ecp_pot%npot(l) 620 n = potential%ecp_pot%nrpot(k, l) 621 alpha = potential%ecp_pot%bpot(k, l) 622 cpot(:) = cpot + potential%ecp_pot%apot(k, l)*rad**(n - 2)*EXP(-alpha*rad**2) 623 END DO 624 DO i = 1, basis%nbas(l) 625 DO j = i, basis%nbas(l) 626 integrals%hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid) 627 integrals%hnl(j, i, l) = integrals%hnl(i, j, l) 628 END DO 629 END DO 630 END DO 631 DEALLOCATE (cpot) 632 END IF 633 634 END IF 635 636 CALL timestop(handle) 637 638 END SUBROUTINE atom_ppint_setup 639 640! ************************************************************************************************** 641!> \brief ... 642!> \param integrals ... 643!> \param basis ... 644!> \param potential ... 645! ************************************************************************************************** 646 SUBROUTINE upfint_setup(integrals, basis, potential) 647 TYPE(atom_integrals), INTENT(INOUT) :: integrals 648 TYPE(atom_basis_type), INTENT(INOUT) :: basis 649 TYPE(atom_potential_type), INTENT(IN) :: potential 650 651 CHARACTER(len=*), PARAMETER :: routineN = 'upfint_setup', routineP = moduleN//':'//routineN 652 653 CHARACTER(len=4) :: ptype 654 INTEGER :: i, j, k1, k2, la, lb, m, n 655 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: spot 656 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: spmat 657 TYPE(atom_basis_type) :: gbasis 658 659 ! get basis representation on UPF grid 660 CALL atom_basis_gridrep(basis, gbasis, potential%upf_pot%r, potential%upf_pot%rab) 661 662 ! local pseudopotential 663 integrals%core = 0._dp 664 CALL numpot_matrix(integrals%core, potential%upf_pot%vlocal, gbasis, 0) 665 666 ptype = ADJUSTL(TRIM(potential%upf_pot%pseudo_type)) 667 IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN 668 ! non local pseudopotential 669 n = MAXVAL(integrals%n(:)) 670 m = potential%upf_pot%number_of_proj 671 ALLOCATE (spmat(n, m)) 672 spmat = 0.0_dp 673 DO i = 1, m 674 la = potential%upf_pot%lbeta(i) 675 DO j = 1, gbasis%nbas(la) 676 spmat(j, i) = integrate_grid(potential%upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid) 677 END DO 678 END DO 679 DO i = 1, m 680 la = potential%upf_pot%lbeta(i) 681 DO j = 1, m 682 lb = potential%upf_pot%lbeta(j) 683 IF (la == lb) THEN 684 DO k1 = 1, gbasis%nbas(la) 685 DO k2 = 1, gbasis%nbas(la) 686 integrals%hnl(k1, k2, la) = integrals%hnl(k1, k2, la) + & 687 spmat(k1, i)*potential%upf_pot%dion(i, j)*spmat(k2, j) 688 END DO 689 END DO 690 END IF 691 END DO 692 END DO 693 DEALLOCATE (spmat) 694 ELSE IF (ptype(1:2) == "SL") THEN 695 ! semi local pseudopotential 696 DO la = 0, potential%upf_pot%l_max 697 IF (la == potential%upf_pot%l_local) CYCLE 698 m = SIZE(potential%upf_pot%vsemi(:, la + 1)) 699 ALLOCATE (spot(m)) 700 spot(:) = potential%upf_pot%vsemi(:, la + 1) - potential%upf_pot%vlocal(:) 701 n = basis%nbas(la) 702 DO i = 1, n 703 DO j = i, n 704 integrals%core(i, j, la) = integrals%core(i, j, la) + & 705 integrate_grid(spot(:), & 706 gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid) 707 integrals%core(j, i, la) = integrals%core(i, j, la) 708 END DO 709 END DO 710 DEALLOCATE (spot) 711 END DO 712 ELSE 713 CPABORT("Pseudopotential type: ["//ADJUSTL(TRIM(ptype))//"] not known") 714 END IF 715 716 ! release basis representation on UPF grid 717 CALL release_atom_basis(gbasis) 718 719 END SUBROUTINE upfint_setup 720 721! ************************************************************************************************** 722!> \brief ... 723!> \param integrals ... 724!> \param basis ... 725!> \param potential ... 726! ************************************************************************************************** 727 SUBROUTINE sgpint_setup(integrals, basis, potential) 728 TYPE(atom_integrals), INTENT(INOUT) :: integrals 729 TYPE(atom_basis_type), INTENT(INOUT) :: basis 730 TYPE(atom_potential_type), INTENT(IN) :: potential 731 732 CHARACTER(len=*), PARAMETER :: routineN = 'sgpint_setup', routineP = moduleN//':'//routineN 733 734 INTEGER :: i, ia, j, l, m, n, na 735 REAL(KIND=dp) :: a, c, rc, zval 736 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, pgauss 737 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qmat 738 REAL(KIND=dp), DIMENSION(:), POINTER :: rad 739 740 rad => basis%grid%rad 741 m = basis%grid%nr 742 743 ! local pseudopotential 744 integrals%core = 0._dp 745 ALLOCATE (cpot(m)) 746 cpot = 0.0_dp 747 zval = potential%sgp_pot%zion 748 DO i = 1, m 749 rc = rad(i)/potential%sgp_pot%ac_local/SQRT(2.0_dp) 750 cpot(i) = cpot(i) - zval/rad(i)*erf(rc) 751 END DO 752 DO i = 1, potential%sgp_pot%n_local 753 cpot(:) = cpot(:) + potential%sgp_pot%c_local(i)*EXP(-potential%sgp_pot%a_local(i)*rad(:)**2) 754 END DO 755 CALL numpot_matrix(integrals%core, cpot, basis, 0) 756 DEALLOCATE (cpot) 757 758 ! nonlocal pseudopotential 759 integrals%hnl = 0.0_dp 760 IF (potential%sgp_pot%has_nonlocal) THEN 761 ALLOCATE (pgauss(1:m)) 762 n = potential%sgp_pot%n_nonlocal 763 ! 764 DO l = 0, potential%sgp_pot%lmax 765 CPASSERT(l <= UBOUND(basis%nbas, 1)) 766 IF (.NOT. potential%sgp_pot%is_nonlocal(l)) CYCLE 767 ! overlap (a|p) 768 na = basis%nbas(l) 769 ALLOCATE (qmat(na, n)) 770 DO i = 1, n 771 pgauss(:) = 0.0_dp 772 DO j = 1, n 773 a = potential%sgp_pot%a_nonlocal(j) 774 c = potential%sgp_pot%c_nonlocal(j, i, l) 775 pgauss(:) = pgauss(:) + c*EXP(-a*rad(:)**2)*rad(:)**l 776 END DO 777 DO ia = 1, na 778 qmat(ia, i) = SUM(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:)) 779 END DO 780 END DO 781 DO i = 1, na 782 DO j = i, na 783 DO ia = 1, n 784 integrals%hnl(i, j, l) = integrals%hnl(i, j, l) & 785 + qmat(i, ia)*qmat(j, ia)*potential%sgp_pot%h_nonlocal(ia, l) 786 END DO 787 integrals%hnl(j, i, l) = integrals%hnl(i, j, l) 788 END DO 789 END DO 790 DEALLOCATE (qmat) 791 END DO 792 DEALLOCATE (pgauss) 793 END IF 794 795 END SUBROUTINE sgpint_setup 796 797! ************************************************************************************************** 798!> \brief ... 799!> \param integrals ... 800!> \param basis ... 801!> \param reltyp ... 802!> \param zcore ... 803!> \param alpha ... 804! ************************************************************************************************** 805 SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha) 806 TYPE(atom_integrals), INTENT(INOUT) :: integrals 807 TYPE(atom_basis_type), INTENT(INOUT) :: basis 808 INTEGER, INTENT(IN) :: reltyp 809 REAL(dp), INTENT(IN) :: zcore 810 REAL(dp), INTENT(IN), OPTIONAL :: alpha 811 812 CHARACTER(len=*), PARAMETER :: routineN = 'atom_relint_setup', & 813 routineP = moduleN//':'//routineN 814 815 INTEGER :: dkhorder, handle, i, k1, k2, l, m, n, nl 816 REAL(dp) :: ascal 817 REAL(dp), ALLOCATABLE, DIMENSION(:) :: cpot 818 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: modpot 819 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ener, sps 820 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hmat, pvp, sp, tp, vp, wfn 821 822 CALL timeset(routineN, handle) 823 824 SELECT CASE (reltyp) 825 CASE DEFAULT 826 CPABORT("") 827 CASE (do_nonrel_atom, do_zoramp_atom, do_sczoramp_atom) 828 dkhorder = -1 829 CASE (do_dkh0_atom) 830 dkhorder = 0 831 CASE (do_dkh1_atom) 832 dkhorder = 1 833 CASE (do_dkh2_atom) 834 dkhorder = 2 835 CASE (do_dkh3_atom) 836 dkhorder = 3 837 END SELECT 838 839 SELECT CASE (reltyp) 840 CASE DEFAULT 841 CPABORT("") 842 CASE (do_nonrel_atom) 843 ! nothing to do 844 NULLIFY (integrals%tzora, integrals%hdkh) 845 CASE (do_zoramp_atom, do_sczoramp_atom) 846 847 NULLIFY (integrals%hdkh) 848 849 IF (integrals%zorastat == 0) THEN 850 n = MAXVAL(basis%nbas) 851 ALLOCATE (integrals%tzora(n, n, 0:lmat)) 852 integrals%tzora = 0._dp 853 m = basis%grid%nr 854 ALLOCATE (modpot(1:m), cpot(1:m)) 855 CALL calculate_model_potential(modpot, basis%grid, zcore) 856 ! Zora potential 857 cpot(1:m) = modpot(1:m)/(4._dp*c_light_au*c_light_au - 2._dp*modpot(1:m)) 858 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m) 859 CALL numpot_matrix(integrals%tzora, cpot, basis, 0) 860 DO l = 0, lmat 861 nl = basis%nbas(l) 862 integrals%tzora(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*integrals%tzora(1:nl, 1:nl, l) 863 END DO 864 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m) 865 CALL numpot_matrix(integrals%tzora, cpot, basis, 2) 866 ! 867 ! scaled ZORA 868 IF (reltyp == do_sczoramp_atom) THEN 869 ALLOCATE (hmat(n, n, 0:lmat), wfn(n, n, 0:lmat), ener(n, 0:lmat), pvp(n, n, 0:lmat), sps(n, n)) 870 hmat(:, :, :) = integrals%kin + integrals%tzora 871 ! model potential 872 CALL numpot_matrix(hmat, modpot, basis, 0) 873 ! eigenvalues and eigenvectors 874 CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne, lmat) 875 ! relativistic kinetic energy 876 cpot(1:m) = c_light_au*c_light_au/(2._dp*c_light_au*c_light_au - modpot(1:m))**2 877 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m) 878 pvp = 0.0_dp 879 CALL numpot_matrix(pvp, cpot, basis, 0) 880 DO l = 0, lmat 881 nl = basis%nbas(l) 882 pvp(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*pvp(1:nl, 1:nl, l) 883 END DO 884 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m) 885 CALL numpot_matrix(pvp, cpot, basis, 2) 886 ! calculate psi*pvp*psi and the scaled orbital energies 887 ! actually, we directly calculate the energy difference 888 DO l = 0, lmat 889 nl = basis%nbas(l) 890 DO i = 1, integrals%nne(l) 891 IF (ener(i, l) < 0._dp) THEN 892 ascal = SUM(wfn(1:nl, i, l)*MATMUL(pvp(1:nl, 1:nl, l), wfn(1:nl, i, l))) 893 ener(i, l) = ener(i, l)*ascal/(1.0_dp + ascal) 894 ELSE 895 ener(i, l) = 0.0_dp 896 END IF 897 END DO 898 END DO 899 ! correction term is calculated as a projector 900 hmat = 0.0_dp 901 DO l = 0, lmat 902 nl = basis%nbas(l) 903 DO i = 1, integrals%nne(l) 904 DO k1 = 1, nl 905 DO k2 = 1, nl 906 hmat(k1, k2, l) = hmat(k1, k2, l) + ener(i, l)*wfn(k1, i, l)*wfn(k2, i, l) 907 END DO 908 END DO 909 END DO 910 ! transform with overlap matrix 911 sps(1:nl, 1:nl) = MATMUL(integrals%ovlp(1:nl, 1:nl, l), & 912 MATMUL(hmat(1:nl, 1:nl, l), integrals%ovlp(1:nl, 1:nl, l))) 913 ! add scaling correction to tzora 914 integrals%tzora(1:nl, 1:nl, l) = integrals%tzora(1:nl, 1:nl, l) - sps(1:nl, 1:nl) 915 END DO 916 917 DEALLOCATE (hmat, wfn, ener, pvp, sps) 918 END IF 919 ! 920 DEALLOCATE (modpot, cpot) 921 922 integrals%zorastat = 1 923 924 END IF 925 926 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom) 927 928 NULLIFY (integrals%tzora) 929 930 IF (integrals%dkhstat == 0) THEN 931 n = MAXVAL(basis%nbas) 932 ALLOCATE (integrals%hdkh(n, n, 0:lmat)) 933 integrals%hdkh = 0._dp 934 935 m = MAXVAL(basis%nprim) 936 ALLOCATE (tp(m, m, 0:lmat), sp(m, m, 0:lmat), vp(m, m, 0:lmat), pvp(m, m, 0:lmat)) 937 tp = 0._dp 938 sp = 0._dp 939 vp = 0._dp 940 pvp = 0._dp 941 942 SELECT CASE (basis%basis_type) 943 CASE DEFAULT 944 CPABORT("") 945 CASE (GTO_BASIS, CGTO_BASIS) 946 947 DO l = 0, lmat 948 m = basis%nprim(l) 949 IF (m > 0) THEN 950 CALL sg_kinetic(tp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l)) 951 CALL sg_overlap(sp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l)) 952 IF (PRESENT(alpha)) THEN 953 CALL sg_erfc(vp(1:m, 1:m, l), l, alpha, basis%am(1:m, l), basis%am(1:m, l)) 954 ELSE 955 CALL sg_nuclear(vp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l)) 956 END IF 957 CALL sg_kinnuc(pvp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l)) 958 vp(1:m, 1:m, l) = -zcore*vp(1:m, 1:m, l) 959 pvp(1:m, 1:m, l) = -zcore*pvp(1:m, 1:m, l) 960 END IF 961 END DO 962 963 CASE (STO_BASIS) 964 CPABORT("") 965 CASE (NUM_BASIS) 966 CPABORT("") 967 END SELECT 968 969 CALL dkh_integrals(integrals, basis, dkhorder, sp, tp, vp, pvp) 970 971 integrals%dkhstat = 1 972 973 DEALLOCATE (tp, sp, vp, pvp) 974 975 ELSE 976 CPASSERT(ASSOCIATED(integrals%hdkh)) 977 END IF 978 979 END SELECT 980 981 CALL timestop(handle) 982 983 END SUBROUTINE atom_relint_setup 984 985! ************************************************************************************************** 986!> \brief ... 987!> \param integrals ... 988!> \param basis ... 989!> \param order ... 990!> \param sp ... 991!> \param tp ... 992!> \param vp ... 993!> \param pvp ... 994! ************************************************************************************************** 995 SUBROUTINE dkh_integrals(integrals, basis, order, sp, tp, vp, pvp) 996 TYPE(atom_integrals), INTENT(INOUT) :: integrals 997 TYPE(atom_basis_type), INTENT(INOUT) :: basis 998 INTEGER, INTENT(IN) :: order 999 REAL(dp), DIMENSION(:, :, 0:) :: sp, tp, vp, pvp 1000 1001 CHARACTER(len=*), PARAMETER :: routineN = 'dkh_integrals', routineP = moduleN//':'//routineN 1002 1003 INTEGER :: l, m, n 1004 REAL(dp), DIMENSION(:, :, :), POINTER :: hdkh 1005 1006 CPASSERT(order >= 0) 1007 1008 hdkh => integrals%hdkh 1009 1010 DO l = 0, lmat 1011 n = integrals%n(l) 1012 m = basis%nprim(l) 1013 IF (n > 0) THEN 1014 CALL dkh_atom_transformation(sp(1:m, 1:m, l), vp(1:m, 1:m, l), tp(1:m, 1:m, l), pvp(1:m, 1:m, l), m, order) 1015 SELECT CASE (basis%basis_type) 1016 CASE DEFAULT 1017 CPABORT("") 1018 CASE (GTO_BASIS) 1019 CPASSERT(n == m) 1020 integrals%hdkh(1:n, 1:n, l) = tp(1:n, 1:n, l) + vp(1:n, 1:n, l) 1021 CASE (CGTO_BASIS) 1022 CALL contract2(integrals%hdkh(1:n, 1:n, l), tp(1:m, 1:m, l), basis%cm(1:m, 1:n, l)) 1023 CALL contract2add(integrals%hdkh(1:n, 1:n, l), vp(1:m, 1:m, l), basis%cm(1:m, 1:n, l)) 1024 CASE (STO_BASIS) 1025 CPABORT("") 1026 CASE (NUM_BASIS) 1027 CPABORT("") 1028 END SELECT 1029 ELSE 1030 integrals%hdkh(1:n, 1:n, l) = 0._dp 1031 END IF 1032 END DO 1033 1034 END SUBROUTINE dkh_integrals 1035 1036! ************************************************************************************************** 1037!> \brief Calculate overlap matrix between two atomic basis sets. 1038!> \param ovlap overlap matrix <chi_{a,l} | chi_{b,l}> 1039!> \param basis_a first basis set 1040!> \param basis_b second basis set 1041! ************************************************************************************************** 1042 SUBROUTINE atom_basis_projection_overlap(ovlap, basis_a, basis_b) 1043 REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: ovlap 1044 TYPE(atom_basis_type), INTENT(IN) :: basis_a, basis_b 1045 1046 INTEGER :: i, j, l, ma, mb, na, nb 1047 LOGICAL :: ebas 1048 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat 1049 1050 ebas = (basis_a%basis_type == basis_b%basis_type) 1051 1052 ovlap = 0.0_dp 1053 1054 IF (ebas) THEN 1055 SELECT CASE (basis_a%basis_type) 1056 CASE DEFAULT 1057 CPABORT("") 1058 CASE (GTO_BASIS) 1059 DO l = 0, lmat 1060 na = basis_a%nbas(l) 1061 nb = basis_b%nbas(l) 1062 CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l)) 1063 END DO 1064 CASE (CGTO_BASIS) 1065 DO l = 0, lmat 1066 na = basis_a%nbas(l) 1067 nb = basis_b%nbas(l) 1068 ma = basis_a%nprim(l) 1069 mb = basis_b%nprim(l) 1070 ALLOCATE (omat(ma, mb)) 1071 CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l)) 1072 ovlap(1:na, 1:nb, l) = MATMUL(TRANSPOSE(basis_a%cm(1:ma, 1:na, l)), & 1073 MATMUL(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l))) 1074 DEALLOCATE (omat) 1075 END DO 1076 CASE (STO_BASIS) 1077 DO l = 0, lmat 1078 na = basis_a%nbas(l) 1079 nb = basis_b%nbas(l) 1080 CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), & 1081 basis_a%ns(1:na, l), basis_b%as(1:nb, l)) 1082 END DO 1083 CASE (NUM_BASIS) 1084 CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid)) 1085 DO l = 0, lmat 1086 na = basis_a%nbas(l) 1087 nb = basis_b%nbas(l) 1088 DO j = 1, nb 1089 DO i = 1, na 1090 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid) 1091 END DO 1092 END DO 1093 END DO 1094 END SELECT 1095 ELSE 1096 CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid)) 1097 DO l = 0, lmat 1098 na = basis_a%nbas(l) 1099 nb = basis_b%nbas(l) 1100 DO j = 1, nb 1101 DO i = 1, na 1102 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid) 1103 END DO 1104 END DO 1105 END DO 1106 END IF 1107 1108 END SUBROUTINE atom_basis_projection_overlap 1109 1110! ************************************************************************************************** 1111!> \brief Release memory allocated for atomic integrals (valence electrons). 1112!> \param integrals atomic integrals 1113! ************************************************************************************************** 1114 SUBROUTINE atom_int_release(integrals) 1115 TYPE(atom_integrals), INTENT(INOUT) :: integrals 1116 1117 CHARACTER(len=*), PARAMETER :: routineN = 'atom_int_release', & 1118 routineP = moduleN//':'//routineN 1119 1120 INTEGER :: ll 1121 1122 IF (ASSOCIATED(integrals%ovlp)) THEN 1123 DEALLOCATE (integrals%ovlp) 1124 END IF 1125 IF (ASSOCIATED(integrals%kin)) THEN 1126 DEALLOCATE (integrals%kin) 1127 END IF 1128 IF (ASSOCIATED(integrals%conf)) THEN 1129 DEALLOCATE (integrals%conf) 1130 END IF 1131 DO ll = 1, SIZE(integrals%ceri) 1132 IF (ASSOCIATED(integrals%ceri(ll)%int)) THEN 1133 DEALLOCATE (integrals%ceri(ll)%int) 1134 END IF 1135 IF (ASSOCIATED(integrals%eeri(ll)%int)) THEN 1136 DEALLOCATE (integrals%eeri(ll)%int) 1137 END IF 1138 END DO 1139 IF (ASSOCIATED(integrals%utrans)) THEN 1140 DEALLOCATE (integrals%utrans) 1141 END IF 1142 IF (ASSOCIATED(integrals%uptrans)) THEN 1143 DEALLOCATE (integrals%uptrans) 1144 END IF 1145 1146 integrals%status = 0 1147 1148 END SUBROUTINE atom_int_release 1149 1150! ************************************************************************************************** 1151!> \brief Release memory allocated for atomic integrals (core electrons). 1152!> \param integrals atomic integrals 1153! ************************************************************************************************** 1154 SUBROUTINE atom_ppint_release(integrals) 1155 TYPE(atom_integrals), INTENT(INOUT) :: integrals 1156 1157 CHARACTER(len=*), PARAMETER :: routineN = 'atom_ppint_release', & 1158 routineP = moduleN//':'//routineN 1159 1160 IF (ASSOCIATED(integrals%hnl)) THEN 1161 DEALLOCATE (integrals%hnl) 1162 END IF 1163 IF (ASSOCIATED(integrals%core)) THEN 1164 DEALLOCATE (integrals%core) 1165 END IF 1166 IF (ASSOCIATED(integrals%clsd)) THEN 1167 DEALLOCATE (integrals%clsd) 1168 END IF 1169 1170 integrals%ppstat = 0 1171 1172 END SUBROUTINE atom_ppint_release 1173 1174! ************************************************************************************************** 1175!> \brief Release memory allocated for atomic integrals (relativistic effects). 1176!> \param integrals atomic integrals 1177! ************************************************************************************************** 1178 SUBROUTINE atom_relint_release(integrals) 1179 TYPE(atom_integrals), INTENT(INOUT) :: integrals 1180 1181 CHARACTER(len=*), PARAMETER :: routineN = 'atom_relint_release', & 1182 routineP = moduleN//':'//routineN 1183 1184 IF (ASSOCIATED(integrals%tzora)) THEN 1185 DEALLOCATE (integrals%tzora) 1186 END IF 1187 IF (ASSOCIATED(integrals%hdkh)) THEN 1188 DEALLOCATE (integrals%hdkh) 1189 END IF 1190 1191 integrals%zorastat = 0 1192 integrals%dkhstat = 0 1193 1194 END SUBROUTINE atom_relint_release 1195 1196! ************************************************************************************************** 1197!> \brief Calculate model potential. It is useful to guess atomic orbitals. 1198!> \param modpot model potential 1199!> \param grid atomic radial grid 1200!> \param zcore nuclear charge 1201! ************************************************************************************************** 1202 SUBROUTINE calculate_model_potential(modpot, grid, zcore) 1203 REAL(dp), DIMENSION(:), INTENT(INOUT) :: modpot 1204 TYPE(grid_atom_type), INTENT(IN) :: grid 1205 REAL(dp), INTENT(IN) :: zcore 1206 1207 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_model_potential', & 1208 routineP = moduleN//':'//routineN 1209 1210 INTEGER :: i, ii, l, ll, n, z 1211 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pot, rho 1212 TYPE(atom_state) :: state 1213 1214 n = SIZE(modpot) 1215 ALLOCATE (rho(n), pot(n)) 1216 1217 ! fill default occupation 1218 state%core = 0._dp 1219 state%occ = 0._dp 1220 state%multiplicity = -1 1221 z = NINT(zcore) 1222 DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1)) 1223 IF (ptable(z)%e_conv(l) /= 0) THEN 1224 state%maxl_occ = l 1225 ll = 2*(2*l + 1) 1226 DO i = 1, SIZE(state%occ, 2) 1227 ii = ptable(z)%e_conv(l) - (i - 1)*ll 1228 IF (ii <= ll) THEN 1229 state%occ(l, i) = ii 1230 EXIT 1231 ELSE 1232 state%occ(l, i) = ll 1233 END IF 1234 END DO 1235 END IF 1236 END DO 1237 1238 modpot = -zcore/grid%rad(:) 1239 1240 ! Coulomb potential 1241 CALL slater_density(rho, pot, NINT(zcore), state, grid) 1242 CALL coulomb_potential_numeric(pot, rho, grid) 1243 modpot = modpot + pot 1244 1245 ! XC potential 1246 CALL wigner_slater_functional(rho, pot) 1247 modpot = modpot + pot 1248 1249 DEALLOCATE (rho, pot) 1250 1251 END SUBROUTINE calculate_model_potential 1252 1253END MODULE atom_operators 1254