38 #include "./base/base_uses.f90"
44 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'atom_operators'
63 eri_coulomb, eri_exchange, all_nu)
64 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
65 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
66 TYPE(atom_potential_type),
INTENT(IN),
OPTIONAL :: potential
67 LOGICAL,
INTENT(IN),
OPTIONAL :: eri_coulomb, eri_exchange, all_nu
69 CHARACTER(len=*),
PARAMETER :: routinen =
'atom_int_setup'
71 INTEGER :: handle, i, ii, info, ipiv(1000), l, l1, &
72 l2, ll, lwork, m, m1, m2, mm1, mm2, n, &
73 n1, n2, nn1, nn2, nu, nx
74 REAL(kind=
dp) :: om, rc, ron, sc, x
75 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cpot, w, work
76 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: omat, vmat
77 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: eri
79 CALL timeset(routinen, handle)
81 IF (integrals%status == 0)
THEN
82 n = maxval(basis%nbas)
83 integrals%n = basis%nbas
85 IF (
PRESENT(eri_coulomb))
THEN
86 integrals%eri_coulomb = eri_coulomb
88 integrals%eri_coulomb = .false.
90 IF (
PRESENT(eri_exchange))
THEN
91 integrals%eri_exchange = eri_exchange
93 integrals%eri_exchange = .false.
95 IF (
PRESENT(all_nu))
THEN
96 integrals%all_nu = all_nu
98 integrals%all_nu = .false.
101 NULLIFY (integrals%ovlp, integrals%kin, integrals%core, integrals%conf)
102 DO ll = 1,
SIZE(integrals%ceri)
103 NULLIFY (integrals%ceri(ll)%int, integrals%eeri(ll)%int)
106 ALLOCATE (integrals%ovlp(n, n, 0:
lmat))
107 integrals%ovlp = 0._dp
109 ALLOCATE (integrals%kin(n, n, 0:
lmat))
110 integrals%kin = 0._dp
114 IF (
PRESENT(potential))
THEN
115 IF (potential%confinement)
THEN
116 ALLOCATE (integrals%conf(n, n, 0:
lmat))
117 integrals%conf = 0._dp
120 IF (potential%conf_type ==
poly_conf)
THEN
123 cpot(1:m) = (basis%grid%rad(1:m)/rc)**sc
129 IF (basis%grid%rad(i) < ron)
THEN
131 ELSEIF (basis%grid%rad(i) < rc)
THEN
132 x = (basis%grid%rad(i) - ron)/om
134 cpot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
135 x = (rc - basis%grid%rad(i))**2/om/(basis%grid%rad(i) - ron)
149 SELECT CASE (basis%basis_type)
155 CALL sg_overlap(integrals%ovlp(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
156 CALL sg_kinetic(integrals%kin(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
158 IF (integrals%eri_coulomb)
THEN
162 nn1 = (n1*(n1 + 1))/2
165 nn2 = (n2*(n2 + 1))/2
166 IF (integrals%all_nu)
THEN
173 cpassert(ll <=
SIZE(integrals%ceri))
174 ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
175 integrals%ceri(ll)%int = 0._dp
176 eri => integrals%ceri(ll)%int
177 CALL sg_coulomb(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
182 IF (integrals%eri_exchange)
THEN
186 nn1 = (n1*(n1 + 1))/2
189 nn2 = (n2*(n2 + 1))/2
190 DO nu = abs(l1 - l2), l1 + l2, 2
192 cpassert(ll <=
SIZE(integrals%eeri))
193 ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
194 integrals%eeri(ll)%int = 0._dp
195 eri => integrals%eeri(ll)%int
196 CALL sg_exchange(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
205 IF (n > 0 .AND. m > 0)
THEN
206 ALLOCATE (omat(m, m))
208 CALL sg_overlap(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
209 CALL contract2(integrals%ovlp(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
210 CALL sg_kinetic(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
211 CALL contract2(integrals%kin(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
215 IF (integrals%eri_coulomb)
THEN
219 nn1 = (n1*(n1 + 1))/2
221 mm1 = (m1*(m1 + 1))/2
224 nn2 = (n2*(n2 + 1))/2
226 mm2 = (m2*(m2 + 1))/2
227 IF (integrals%all_nu)
THEN
234 cpassert(ll <=
SIZE(integrals%ceri))
235 ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
236 integrals%ceri(ll)%int = 0._dp
237 ALLOCATE (omat(mm1, mm2))
239 eri => integrals%ceri(ll)%int
240 CALL sg_coulomb(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
241 CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))
248 IF (integrals%eri_exchange)
THEN
252 nn1 = (n1*(n1 + 1))/2
254 mm1 = (m1*(m1 + 1))/2
257 nn2 = (n2*(n2 + 1))/2
259 mm2 = (m2*(m2 + 1))/2
260 DO nu = abs(l1 - l2), l1 + l2, 2
262 cpassert(ll <=
SIZE(integrals%eeri))
263 ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
264 integrals%eeri(ll)%int = 0._dp
265 ALLOCATE (omat(mm1, mm2))
267 eri => integrals%eeri(ll)%int
268 CALL sg_exchange(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
269 CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))
279 CALL sto_overlap(integrals%ovlp(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
280 basis%ns(1:n, l), basis%as(1:n, l))
281 CALL sto_kinetic(integrals%kin(1:n, 1:n, l), l, basis%ns(1:n, l), basis%as(1:n, l), &
282 basis%ns(1:n, l), basis%as(1:n, l))
284 cpassert(.NOT. integrals%eri_coulomb)
285 cpassert(.NOT. integrals%eri_exchange)
291 NULLIFY (integrals%utrans, integrals%uptrans)
292 n = maxval(basis%nbas)
293 ALLOCATE (integrals%utrans(n, n, 0:
lmat), integrals%uptrans(n, n, 0:
lmat))
294 integrals%utrans = 0._dp
295 integrals%uptrans = 0._dp
296 integrals%nne = integrals%n
298 ALLOCATE (omat(n, n), vmat(n, n), w(n), work(lwork))
302 omat(1:n, 1:n) = integrals%ovlp(1:n, 1:n, l)
303 CALL jacobi(omat(1:n, 1:n), w(1:n), vmat(1:n, 1:n))
304 omat(1:n, 1:n) = vmat(1:n, 1:n)
307 IF (w(i) > basis%eps_eig)
THEN
309 integrals%utrans(1:n, ii, l) = omat(1:n, i)/sqrt(w(i))
312 integrals%nne(l) = ii
314 omat(1:ii, 1:ii) = matmul(transpose(integrals%utrans(1:n, 1:ii, l)), integrals%utrans(1:n, 1:ii, l))
316 integrals%uptrans(i, i, l) = 1._dp
318 CALL lapack_sgesv(ii, ii, omat(1:ii, 1:ii), ii, ipiv, integrals%uptrans(1:ii, 1:ii, l), ii, info)
323 DEALLOCATE (omat, vmat, w, work)
326 CALL timestop(handle)
337 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
338 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
339 TYPE(atom_potential_type),
INTENT(IN) :: potential
341 CHARACTER(len=*),
PARAMETER :: routinen =
'atom_ppint_setup'
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
349 CALL timeset(routinen, handle)
351 IF (integrals%ppstat == 0)
THEN
352 n = maxval(basis%nbas)
353 integrals%n = basis%nbas
355 NULLIFY (integrals%core, integrals%hnl)
357 ALLOCATE (integrals%hnl(n, n, 0:
lmat))
358 integrals%hnl = 0._dp
360 ALLOCATE (integrals%core(n, n, 0:
lmat))
361 integrals%core = 0._dp
363 ALLOCATE (integrals%clsd(n, n, 0:
lmat))
364 integrals%clsd = 0._dp
368 SELECT CASE (basis%basis_type)
373 SELECT CASE (potential%ppot_type)
374 CASE (no_pseudo, ecp_pseudo)
377 CALL sg_nuclear(integrals%core(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
380 alpha = 1._dp/potential%gth_pot%rc/sqrt(2._dp)
383 ALLOCATE (omat(n, n), spmat(n, 5))
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
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)
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)
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)
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)
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)
418 m = potential%gth_pot%nl(l)
420 CALL sg_proj_ol(spmat(1:n, i), l, basis%am(1:n, l), i - 1, potential%gth_pot%rcnl(l))
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))))
425 DEALLOCATE (omat, spmat)
428 CALL upfint_setup(integrals, basis, potential)
430 CALL sgpint_setup(integrals, basis, potential)
437 SELECT CASE (potential%ppot_type)
438 CASE (no_pseudo, ecp_pseudo)
442 ALLOCATE (omat(m, m))
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))
450 alpha = 1._dp/potential%gth_pot%rc/sqrt(2._dp)
454 IF (n > 0 .AND. m > 0)
THEN
455 ALLOCATE (omat(m, m), spmat(n, 5), xmat(m))
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
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))
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)
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))
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)
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))
491 k = potential%gth_pot%nl(l)
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))
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))))
501 DEALLOCATE (omat, spmat, xmat)
505 CALL upfint_setup(integrals, basis, potential)
507 CALL sgpint_setup(integrals, basis, potential)
514 SELECT CASE (potential%ppot_type)
515 CASE (no_pseudo, ecp_pseudo)
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))
522 rad => basis%grid%rad
525 rc = potential%gth_pot%rc
526 alpha = 1._dp/rc/sqrt(2._dp)
529 integrals%core = 0._dp
531 cpot(i) = potential%gth_pot%zion*erfc(alpha*rad(i))/rad(i)
533 DO i = 1, potential%gth_pot%ncl
535 cpot(1:m) = cpot(1:m) + potential%gth_pot%cl(i)*(rad/rc)**ii*exp(-0.5_dp*(rad/rc)**2)
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)
542 cpot(1:m) = cpot(1:m) + potential%gth_pot%cval_lpot(i, k)*(rad/al)**ii*exp(-0.5_dp*(rad/al)**2)
549 ALLOCATE (omat(n, n))
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)
557 IF (potential%gth_pot%lsdpot)
THEN
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)
563 cpot(:) = cpot + potential%gth_pot%cval_lsd(i, k)*(rad/al)**ii*exp(-0.5_dp*(rad/al)**2)
572 ALLOCATE (spmat(n, 5))
574 k = potential%gth_pot%nl(l)
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)
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))))
590 CALL upfint_setup(integrals, basis, potential)
592 CALL sgpint_setup(integrals, basis, potential)
602 IF (potential%ppot_type == ecp_pseudo)
THEN
604 integrals%core = -potential%ecp_pot%zion*integrals%core
607 rad => basis%grid%rad
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)
617 DO l = 0, min(potential%ecp_pot%lmax,
lmat)
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)
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)
636 CALL timestop(handle)
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
651 CHARACTER(len=4) :: ptype
652 INTEGER :: i, j, k1, k2, la, lb, m, n
653 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: spot
654 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: spmat
655 TYPE(atom_basis_type) :: gbasis
661 integrals%core = 0._dp
662 CALL numpot_matrix(integrals%core, potential%upf_pot%vlocal, gbasis, 0)
664 ptype = adjustl(trim(potential%upf_pot%pseudo_type))
665 IF (ptype(1:2) ==
"NC" .OR. ptype(1:2) ==
"US")
THEN
667 n = maxval(integrals%n(:))
668 m = potential%upf_pot%number_of_proj
669 ALLOCATE (spmat(n, m))
672 la = potential%upf_pot%lbeta(i)
673 DO j = 1, gbasis%nbas(la)
674 spmat(j, i) = integrate_grid(potential%upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
678 la = potential%upf_pot%lbeta(i)
680 lb = potential%upf_pot%lbeta(j)
682 DO k1 = 1, gbasis%nbas(la)
683 DO k2 = 1, gbasis%nbas(la)
684 integrals%hnl(k1, k2, la) = integrals%hnl(k1, k2, la) + &
685 spmat(k1, i)*potential%upf_pot%dion(i, j)*spmat(k2, j)
692 ELSE IF (ptype(1:2) ==
"SL")
THEN
694 DO la = 0, potential%upf_pot%l_max
695 IF (la == potential%upf_pot%l_local) cycle
696 m =
SIZE(potential%upf_pot%vsemi(:, la + 1))
698 spot(:) = potential%upf_pot%vsemi(:, la + 1) - potential%upf_pot%vlocal(:)
702 integrals%core(i, j, la) = integrals%core(i, j, la) + &
703 integrate_grid(spot(:), &
704 gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
705 integrals%core(j, i, la) = integrals%core(i, j, la)
711 cpabort(
"Pseudopotential type: ["//adjustl(trim(ptype))//
"] not known")
717 END SUBROUTINE upfint_setup
725 SUBROUTINE sgpint_setup(integrals, basis, potential)
726 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
727 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
728 TYPE(atom_potential_type),
INTENT(IN) :: potential
730 INTEGER :: i, ia, j, l, m, n, na
731 REAL(kind=
dp) :: a, c, rc, zval
732 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cpot, pgauss
733 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: qmat
734 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rad
736 rad => basis%grid%rad
740 integrals%core = 0._dp
743 zval = potential%sgp_pot%zion
745 rc = rad(i)/potential%sgp_pot%ac_local/sqrt(2.0_dp)
746 cpot(i) = cpot(i) - zval/rad(i)*erf(rc)
748 DO i = 1, potential%sgp_pot%n_local
749 cpot(:) = cpot(:) + potential%sgp_pot%c_local(i)*exp(-potential%sgp_pot%a_local(i)*rad(:)**2)
755 integrals%hnl = 0.0_dp
756 IF (potential%sgp_pot%has_nonlocal)
THEN
757 ALLOCATE (pgauss(1:m))
758 n = potential%sgp_pot%n_nonlocal
760 DO l = 0, potential%sgp_pot%lmax
761 cpassert(l <= ubound(basis%nbas, 1))
762 IF (.NOT. potential%sgp_pot%is_nonlocal(l)) cycle
765 ALLOCATE (qmat(na, n))
769 a = potential%sgp_pot%a_nonlocal(j)
770 c = potential%sgp_pot%c_nonlocal(j, i, l)
771 pgauss(:) = pgauss(:) + c*exp(-a*rad(:)**2)*rad(:)**l
774 qmat(ia, i) = sum(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:))
780 integrals%hnl(i, j, l) = integrals%hnl(i, j, l) &
781 + qmat(i, ia)*qmat(j, ia)*potential%sgp_pot%h_nonlocal(ia, l)
783 integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
791 END SUBROUTINE sgpint_setup
802 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
803 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
804 INTEGER,
INTENT(IN) :: reltyp
805 REAL(
dp),
INTENT(IN) :: zcore
806 REAL(
dp),
INTENT(IN),
OPTIONAL :: alpha
808 CHARACTER(len=*),
PARAMETER :: routinen =
'atom_relint_setup'
810 INTEGER :: dkhorder, handle, i, k1, k2, l, m, n, nl
812 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: cpot
813 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: modpot
814 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ener, sps
815 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: hmat, pvp, sp, tp, vp, wfn
817 CALL timeset(routinen, handle)
839 NULLIFY (integrals%tzora, integrals%hdkh)
842 NULLIFY (integrals%hdkh)
844 IF (integrals%zorastat == 0)
THEN
845 n = maxval(basis%nbas)
846 ALLOCATE (integrals%tzora(n, n, 0:
lmat))
847 integrals%tzora = 0._dp
849 ALLOCATE (modpot(1:m), cpot(1:m))
853 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
857 integrals%tzora(1:nl, 1:nl, l) = real(l*(l + 1),
dp)*integrals%tzora(1:nl, 1:nl, l)
859 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
864 ALLOCATE (hmat(n, n, 0:
lmat), wfn(n, n, 0:
lmat), ener(n, 0:
lmat), pvp(n, n, 0:
lmat), sps(n, n))
865 hmat(:, :, :) = integrals%kin + integrals%tzora
869 CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne,
lmat)
872 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
877 pvp(1:nl, 1:nl, l) = real(l*(l + 1),
dp)*pvp(1:nl, 1:nl, l)
879 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
885 DO i = 1, integrals%nne(l)
886 IF (ener(i, l) < 0._dp)
THEN
887 ascal = sum(wfn(1:nl, i, l)*matmul(pvp(1:nl, 1:nl, l), wfn(1:nl, i, l)))
888 ener(i, l) = ener(i, l)*ascal/(1.0_dp + ascal)
898 DO i = 1, integrals%nne(l)
901 hmat(k1, k2, l) = hmat(k1, k2, l) + ener(i, l)*wfn(k1, i, l)*wfn(k2, i, l)
906 sps(1:nl, 1:nl) = matmul(integrals%ovlp(1:nl, 1:nl, l), &
907 matmul(hmat(1:nl, 1:nl, l), integrals%ovlp(1:nl, 1:nl, l)))
909 integrals%tzora(1:nl, 1:nl, l) = integrals%tzora(1:nl, 1:nl, l) - sps(1:nl, 1:nl)
912 DEALLOCATE (hmat, wfn, ener, pvp, sps)
915 DEALLOCATE (modpot, cpot)
917 integrals%zorastat = 1
923 NULLIFY (integrals%tzora)
925 IF (integrals%dkhstat == 0)
THEN
926 n = maxval(basis%nbas)
927 ALLOCATE (integrals%hdkh(n, n, 0:
lmat))
928 integrals%hdkh = 0._dp
930 m = maxval(basis%nprim)
931 ALLOCATE (tp(m, m, 0:
lmat), sp(m, m, 0:
lmat), vp(m, m, 0:
lmat), pvp(m, m, 0:
lmat))
937 SELECT CASE (basis%basis_type)
945 CALL sg_kinetic(tp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
946 CALL sg_overlap(sp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
947 IF (
PRESENT(alpha))
THEN
948 CALL sg_erfc(vp(1:m, 1:m, l), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
950 CALL sg_nuclear(vp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
952 CALL sg_kinnuc(pvp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
953 vp(1:m, 1:m, l) = -zcore*vp(1:m, 1:m, l)
954 pvp(1:m, 1:m, l) = -zcore*pvp(1:m, 1:m, l)
964 CALL dkh_integrals(integrals, basis, dkhorder, sp, tp, vp, pvp)
966 integrals%dkhstat = 1
968 DEALLOCATE (tp, sp, vp, pvp)
971 cpassert(
ASSOCIATED(integrals%hdkh))
976 CALL timestop(handle)
990 SUBROUTINE dkh_integrals(integrals, basis, order, sp, tp, vp, pvp)
991 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
992 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
993 INTEGER,
INTENT(IN) :: order
994 REAL(
dp),
DIMENSION(:, :, 0:) :: sp, tp, vp, pvp
997 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: hdkh
1001 hdkh => integrals%hdkh
1007 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)
1008 SELECT CASE (basis%basis_type)
1013 integrals%hdkh(1:n, 1:n, l) = tp(1:n, 1:n, l) + vp(1:n, 1:n, l)
1015 CALL contract2(integrals%hdkh(1:n, 1:n, l), tp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1016 CALL contract2add(integrals%hdkh(1:n, 1:n, l), vp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1023 integrals%hdkh(1:n, 1:n, l) = 0._dp
1027 END SUBROUTINE dkh_integrals
1036 REAL(kind=
dp),
DIMENSION(:, :, 0:),
INTENT(OUT) :: ovlap
1037 TYPE(atom_basis_type),
INTENT(IN) :: basis_a, basis_b
1039 INTEGER :: i, j, l, ma, mb, na, nb
1041 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: omat
1043 ebas = (basis_a%basis_type == basis_b%basis_type)
1048 SELECT CASE (basis_a%basis_type)
1053 na = basis_a%nbas(l)
1054 nb = basis_b%nbas(l)
1055 CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l))
1059 na = basis_a%nbas(l)
1060 nb = basis_b%nbas(l)
1061 ma = basis_a%nprim(l)
1062 mb = basis_b%nprim(l)
1063 ALLOCATE (omat(ma, mb))
1064 CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l))
1065 ovlap(1:na, 1:nb, l) = matmul(transpose(basis_a%cm(1:ma, 1:na, l)), &
1066 matmul(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l)))
1071 na = basis_a%nbas(l)
1072 nb = basis_b%nbas(l)
1073 CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), &
1074 basis_a%ns(1:na, l), basis_b%as(1:nb, l))
1079 na = basis_a%nbas(l)
1080 nb = basis_b%nbas(l)
1083 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1091 na = basis_a%nbas(l)
1092 nb = basis_b%nbas(l)
1095 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1108 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
1112 IF (
ASSOCIATED(integrals%ovlp))
THEN
1113 DEALLOCATE (integrals%ovlp)
1115 IF (
ASSOCIATED(integrals%kin))
THEN
1116 DEALLOCATE (integrals%kin)
1118 IF (
ASSOCIATED(integrals%conf))
THEN
1119 DEALLOCATE (integrals%conf)
1121 DO ll = 1,
SIZE(integrals%ceri)
1122 IF (
ASSOCIATED(integrals%ceri(ll)%int))
THEN
1123 DEALLOCATE (integrals%ceri(ll)%int)
1125 IF (
ASSOCIATED(integrals%eeri(ll)%int))
THEN
1126 DEALLOCATE (integrals%eeri(ll)%int)
1129 IF (
ASSOCIATED(integrals%utrans))
THEN
1130 DEALLOCATE (integrals%utrans)
1132 IF (
ASSOCIATED(integrals%uptrans))
THEN
1133 DEALLOCATE (integrals%uptrans)
1136 integrals%status = 0
1145 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
1147 IF (
ASSOCIATED(integrals%hnl))
THEN
1148 DEALLOCATE (integrals%hnl)
1150 IF (
ASSOCIATED(integrals%core))
THEN
1151 DEALLOCATE (integrals%core)
1153 IF (
ASSOCIATED(integrals%clsd))
THEN
1154 DEALLOCATE (integrals%clsd)
1157 integrals%ppstat = 0
1166 TYPE(atom_integrals),
INTENT(INOUT) :: integrals
1168 IF (
ASSOCIATED(integrals%tzora))
THEN
1169 DEALLOCATE (integrals%tzora)
1171 IF (
ASSOCIATED(integrals%hdkh))
THEN
1172 DEALLOCATE (integrals%hdkh)
1175 integrals%zorastat = 0
1176 integrals%dkhstat = 0
1187 REAL(
dp),
DIMENSION(:),
INTENT(INOUT) :: modpot
1188 TYPE(grid_atom_type),
INTENT(IN) :: grid
1189 REAL(
dp),
INTENT(IN) :: zcore
1191 INTEGER :: i, ii, l, ll, n, z
1192 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: pot, rho
1193 TYPE(atom_state) :: state
1196 ALLOCATE (rho(n), pot(n))
1201 state%multiplicity = -1
1203 DO l = 0, min(
lmat, ubound(
ptable(z)%e_conv, 1))
1204 IF (
ptable(z)%e_conv(l) /= 0)
THEN
1207 DO i = 1,
SIZE(state%occ, 2)
1208 ii =
ptable(z)%e_conv(l) - (i - 1)*ll
1210 state%occ(l, i) = ii
1213 state%occ(l, i) = ll
1219 modpot = -zcore/grid%rad(:)
1224 modpot = modpot + pot
1228 modpot = modpot + pot
1230 DEALLOCATE (rho, pot)
subroutine, public sg_kinnuc(umat, l, pa, pb)
...
subroutine, public sg_nuclear(umat, l, pa, pb)
...
subroutine, public sg_coulomb(eri, nu, pa, lab, pc, lcd)
...
subroutine, public sg_exchange(eri, nu, pa, lac, pb, lbd)
...
subroutine, public sg_erfc(umat, l, a, pa, pb)
...
subroutine, public sg_kinetic(kmat, l, pa, pb)
...
subroutine, public sto_nuclear(umat, na, pa, nb, pb)
...
subroutine, public sto_overlap(smat, na, pa, nb, pb)
...
subroutine, public sto_kinetic(kmat, l, na, pa, nb, pb)
...
subroutine, public sg_gpot(vpmat, k, rc, l, pa, pb)
...
subroutine, public sg_proj_ol(spmat, l, p, k, rc)
...
subroutine, public sg_erf(upmat, l, a, pa, pb)
...
subroutine, public sg_overlap(smat, l, pa, pb)
...
Calculate the atomic operator matrices.
subroutine, public atom_ppint_release(integrals)
Release memory allocated for atomic integrals (core electrons).
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
subroutine, public atom_basis_projection_overlap(ovlap, basis_a, basis_b)
Calculate overlap matrix between two atomic basis sets.
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
subroutine, public calculate_model_potential(modpot, grid, zcore)
Calculate model potential. It is useful to guess atomic orbitals.
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Define the atom type and its sub types.
integer, parameter, public num_basis
integer, parameter, public cgto_basis
integer, parameter, public gto_basis
integer, parameter, public sto_basis
logical function, public atom_compare_grids(grid1, grid2)
...
integer, parameter, public lmat
subroutine, public release_atom_basis(basis)
...
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
Some basic routines for atomic calculations.
subroutine, public slater_density(density1, density2, zcore, state, grid)
Calculate Slater density on a radial grid.
subroutine, public contract2(int, omat, cm)
Transform a matrix expressed in terms of a uncontracted basis set to a contracted one.
pure subroutine, public wigner_slater_functional(rho, vxc)
Calculate the functional derivative of the Wigner (correlation) - Slater (exchange) density functiona...
subroutine, public contract4(eri, omat, cm1, cm2)
Contract a matrix of Electron Repulsion Integrals (ERI-s).
subroutine, public numpot_matrix(imat, cpot, basis, derivatives)
Calculate a potential matrix by integrating the potential on an atomic radial grid.
subroutine, public atom_solve(hmat, umat, orb, ener, nb, nv, maxl)
Solve the generalised eigenproblem for every angular momentum.
subroutine, public contract2add(int, omat, cm)
Same as contract2(), but add the new contracted matrix to the old one instead of overwriting it.
subroutine, public coulomb_potential_numeric(cpot, density, grid)
Numerically compute the Coulomb potential on an atomic radial grid.
subroutine, public dkh_atom_transformation(s, v, h, pVp, n, dkh_order)
...
Defines the basic variable types.
integer, parameter, public dp
Interface to the LAPACK F77 library.
Definition of mathematical constants and functions.
real(kind=dp), dimension(0:maxfac), parameter, public gamma1
real(kind=dp), parameter, public sqrt2
Collection of simple mathematical functions and subroutines.
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
real(kind=dp), parameter, public c_light_au