57#include "./base/base_uses.f90"
63 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'atom_kind_orbitals'
92 density, wavefunction, wfninfo, confine, xc_section, nocc, &
93 which_l, which_n, proj_shell_charge, ao_coef)
97 INTEGER,
INTENT(IN),
OPTIONAL :: iunit
98 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
100 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: density
101 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: wavefunction, wfninfo
102 LOGICAL,
INTENT(IN),
OPTIONAL :: confine
104 INTEGER,
DIMENSION(:),
OPTIONAL :: nocc
105 INTEGER,
INTENT(IN),
OPTIONAL :: which_l, which_n
106 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: proj_shell_charge(4)
107 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT), &
108 OPTIONAL,
POINTER :: ao_coef
110 INTEGER :: i, ii, j, k, k1, k2, l, ll, m, mb, mo, &
111 nr, nset, nsgf, projector_index, &
113 INTEGER,
DIMENSION(0:lmat) :: nbb
114 INTEGER,
DIMENSION(0:lmat, 10) :: ncalc, ncore,
nelem
115 INTEGER,
DIMENSION(0:lmat, 100) :: set_index, shell_index
116 INTEGER,
DIMENSION(:),
POINTER :: nshell
117 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgf, ls
118 LOGICAL :: ecp_semi_local, ghost, has_pp, uks
119 REAL(kind=
dp) :: ok, qleft, scal, take, zeff
120 REAL(kind=
dp),
DIMENSION(0:lmat, 10) :: occupation_neutral
121 REAL(kind=
dp),
DIMENSION(0:lmat, 10, 2) :: edelta
132 IF (
PRESENT(ao_coef) .OR.
PRESENT(which_l) .OR. &
133 PRESENT(which_n) .OR.
PRESENT(proj_shell_charge))
THEN
134 IF (.NOT.
PRESENT(ao_coef) .OR. .NOT.
PRESENT(which_l) .OR. &
135 .NOT.
PRESENT(which_n) .OR. .NOT.
PRESENT(proj_shell_charge))
THEN
136 cpabort(
"ao_coef, which_l, which_n, and proj_shell_charge must be provided together")
143 IF (
PRESENT(xc_section))
THEN
144 atom%xc_section => xc_section
146 NULLIFY (
atom%xc_section)
150 NULLIFY (all_potential, gth_potential, sgp_potential, orb_basis_set)
152 basis_set=orb_basis_set, &
154 all_potential=all_potential, &
155 gth_potential=gth_potential, &
156 sgp_potential=sgp_potential)
158 has_pp =
ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential)
170 ALLOCATE (potential, integrals)
172 IF (
PRESENT(confine))
THEN
173 potential%confinement = confine
175 IF (
ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential))
THEN
176 potential%confinement = .true.
178 potential%confinement = .false.
182 potential%acon = 0.1_dp
183 potential%rcon = 2.0_dp*
ptable(z)%vdw_radius*
bohr
184 potential%scon = 2.0_dp
186 IF (
ASSOCIATED(gth_potential))
THEN
190 CALL set_atom(
atom, zcore=nint(zeff), potential=potential)
191 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
192 CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
193 IF (ecp_semi_local)
THEN
195 CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
196 potential%ecp_pot%symbol =
ptable(z)%symbol
199 CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
200 potential%sgp_pot%symbol =
ptable(z)%symbol
203 CALL set_atom(
atom, zcore=nint(zeff), potential=potential)
217 atom%optimization%damping = 0.2_dp
218 atom%optimization%eps_scf = 1.e-6_dp
219 atom%optimization%eps_diis = 100._dp
220 atom%optimization%max_iter = 50
221 atom%optimization%n_diis = 5
230 rks=
PRESENT(ao_coef))
233 IF (sum(abs(edelta)) > 0.0_dp)
THEN
241 ALLOCATE (
atom%state)
243 atom%state%core = 0._dp
245 atom%state%occ = 0._dp
248 edelta(0:
lmat, 1:7, 1) + edelta(0:
lmat, 1:7, 2)
252 atom%state%occupation = 0._dp
256 IF (ncalc(l, i) > 0)
THEN
259 atom%state%occupation(l, k) = real(ncalc(l, i),
dp) + &
260 edelta(l, i, 1) + edelta(l, i, 2)
261 atom%state%occa(l, k) = 0.5_dp*real(ncalc(l, i),
dp) + edelta(l, i, 1)
262 atom%state%occb(l, k) = 0.5_dp*real(ncalc(l, i),
dp) + edelta(l, i, 2)
264 atom%state%occupation(l, k) = real(ncalc(l, i),
dp)
268 ok = real(2*l + 1, kind=
dp)
271 atom%state%occ(l, i) = min(
atom%state%occ(l, i), 2.0_dp*ok)
272 atom%state%occa(l, i) = min(
atom%state%occa(l, i), ok)
273 atom%state%occb(l, i) = min(
atom%state%occb(l, i), ok)
274 atom%state%occupation(l, i) =
atom%state%occa(l, i) +
atom%state%occb(l, i)
278 atom%state%occ(l, i) = min(
atom%state%occ(l, i), 2.0_dp*ok)
279 atom%state%occupation(l, i) = min(
atom%state%occupation(l, i), 2.0_dp*ok)
284 IF (
PRESENT(proj_shell_charge))
THEN
286 occupation_neutral(0:
lmat, 1:10) =
atom%state%occupation(0:
lmat, 1:10)
288 CALL cp_abort(__location__,
"proj_shell_charge in"// &
289 " calculate_atomic_orbitals is implemented only for RKS")
291 DO l = 0, min(3,
lmat)
292 qleft = proj_shell_charge(l + 1)
293 IF (qleft > 0.0_dp)
THEN
296 IF (
atom%state%occupation(l, k) > 0.0_dp)
THEN
297 take = min(qleft,
atom%state%occupation(l, k))
298 atom%state%occupation(l, k) =
atom%state%occupation(l, k) - take
300 IF (qleft <= 1.0e-12_dp)
EXIT
303 IF (qleft > 1.0e-10_dp)
THEN
304 CALL cp_abort(__location__,
"proj_shell_charge removes"// &
305 " more electrons than available in this angular channel")
308 qleft = proj_shell_charge(l + 1)
310 IF (
atom%state%occ(l, k) > 0.0_dp)
THEN
311 take = min(qleft,
atom%state%occ(l, k))
312 atom%state%occ(l, k) =
atom%state%occ(l, k) - take
314 IF (qleft <= 1.0e-12_dp)
EXIT
317 ELSE IF (qleft < 0.0_dp)
THEN
318 CALL cp_abort(__location__,
"Negative proj_shell_charge"// &
319 " is not supported in this implementation")
325 atom%state%multiplicity = nint(abs(sum(
atom%state%occa -
atom%state%occb)) + 1)
327 atom%state%multiplicity = -1
332 IF (
PRESENT(proj_shell_charge))
THEN
336 atom%state%maxl_calc =
atom%state%maxl_occ
337 atom%state%maxn_calc =
atom%state%maxn_occ
341 IF (
PRESENT(nocc) .AND. ghost)
THEN
343 ELSE IF (
PRESENT(nocc))
THEN
348 IF (
atom%state%occa(l, k) > 0.0_dp)
THEN
349 nocc(1) = nocc(1) + 2*l + 1
351 IF (
atom%state%occb(l, k) > 0.0_dp)
THEN
352 nocc(2) = nocc(2) + 2*l + 1
355 IF (
atom%state%occupation(l, k) > 0.0_dp)
THEN
356 nocc(1) = nocc(1) + 2*l + 1
357 nocc(2) = nocc(2) + 2*l + 1
372 NULLIFY (integrals%tzora, integrals%hdkh)
377 mo = maxval(
atom%state%maxn_calc)
378 mb = maxval(
atom%basis%nbas)
382 IF (.NOT. ghost)
THEN
383 IF (
PRESENT(iunit))
THEN
389 IF (
PRESENT(pmat))
THEN
392 nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
404 shell_index(l, k) = j
409 IF (
ASSOCIATED(pmat))
THEN
412 ALLOCATE (pmat(nsgf, nsgf, 2))
414 IF (.NOT. ghost)
THEN
417 DO k1 = 1,
atom%basis%nbas(l)
418 DO k2 = 1,
atom%basis%nbas(l)
419 scal = sqrt(
atom%integrals%ovlp(k1, k1, l)*
atom%integrals%ovlp(k2, k2, l))/real(2*l + 1, kind=
dp)
420 i = first_sgf(shell_index(l, k1), set_index(l, k1))
421 j = first_sgf(shell_index(l, k2), set_index(l, k2))
424 pmat(i + m, j + m, 1) =
atom%orbitals%pmata(k1, k2, l)*scal
425 pmat(i + m, j + m, 2) =
atom%orbitals%pmatb(k1, k2, l)*scal
429 pmat(i + m, j + m, 1) =
atom%orbitals%pmat(k1, k2, l)*scal
436 pmat(:, :, 1) = pmat(:, :, 1) + pmat(:, :, 2)
437 pmat(:, :, 2) = pmat(:, :, 1) - 2.0_dp*pmat(:, :, 2)
442 IF (
PRESENT(fmat))
THEN
446 nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
458 shell_index(l, k) = j
462 IF (uks) cpabort(
"calculate_atomic_orbitals: only RKS is implemented")
463 IF (
ASSOCIATED(fmat)) cpabort(
"fmat already associated")
464 IF (.NOT.
ASSOCIATED(
atom%fmat)) cpabort(
"atom%fmat not associated")
465 ALLOCATE (fmat(nsgf, nsgf, 1))
467 IF (.NOT. ghost)
THEN
470 DO k1 = 1,
atom%basis%nbas(l)
471 DO k2 = 1,
atom%basis%nbas(l)
472 scal = sqrt(
atom%integrals%ovlp(k1, k1, l)*
atom%integrals%ovlp(k2, k2, l))
473 i = first_sgf(shell_index(l, k1), set_index(l, k1))
474 j = first_sgf(shell_index(l, k2), set_index(l, k2))
476 fmat(i + m, j + m, 1) =
atom%fmat%op(k1, k2, l)/scal
486 IF (
PRESENT(density))
THEN
487 IF (
ASSOCIATED(density))
DEALLOCATE (density)
488 ALLOCATE (density(nr))
496 IF (
PRESENT(wavefunction))
THEN
497 cpassert(
PRESENT(wfninfo))
498 IF (
ASSOCIATED(wavefunction))
DEALLOCATE (wavefunction)
499 IF (
ASSOCIATED(wfninfo))
DEALLOCATE (wfninfo)
500 mo = sum(
atom%state%maxn_occ)
501 ALLOCATE (wavefunction(nr, mo), wfninfo(2, mo))
502 wavefunction = 0.0_dp
503 IF (.NOT. ghost)
THEN
506 DO i = 1,
atom%state%maxn_occ(l)
507 IF (
atom%state%occupation(l, i) > 0.0_dp)
THEN
509 wfninfo(1, ii) =
atom%state%occupation(l, i)
510 wfninfo(2, ii) = real(l,
dp)
511 DO j = 1,
atom%basis%nbas(l)
512 wavefunction(:, ii) = wavefunction(:, ii) + &
513 atom%orbitals%wfn(j, i, l)*basis%bf(:, j, l)
522 IF (
PRESENT(ao_coef))
THEN
523 cpassert(
PRESENT(which_l))
524 cpassert(
PRESENT(which_n))
525 cpassert(which_l >= 0 .AND. which_l <=
lmat)
526 valence_n = count(ncore(which_l, :) > 0) + which_l + 1
527 projector_index = which_n - valence_n + 1
528 IF (projector_index < 1 .OR. projector_index >
SIZE(
atom%orbitals%wfn, 2))
THEN
529 CALL cp_abort(__location__, &
530 "The requested tensorial atomic projector N is outside "// &
531 "the available orbital range for the selected L.")
533 IF (
ASSOCIATED(ao_coef))
DEALLOCATE (ao_coef)
534 ALLOCATE (ao_coef(
SIZE(
atom%orbitals%wfn(:, 1, which_l))))
535 ao_coef(:) =
atom%orbitals%wfn(:, projector_index, which_l)
536 IF (maxval(abs(ao_coef)) <= 100.0_dp*epsilon(1.0_dp))
THEN
537 CALL cp_abort(__location__, &
538 "The requested atomic projector is zero. Check the "// &
539 "specified N and L quantum numbers.")
541 DO j = 1,
SIZE(ao_coef(:))
542 ao_coef(j) = ao_coef(j)*sqrt(
atom%integrals%ovlp(j, j, which_l))
554 DEALLOCATE (potential, basis, integrals)
570 optbasis, allelectron, confine)
571 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: density
574 INTEGER,
INTENT(IN) :: ngto
575 INTEGER,
INTENT(IN),
OPTIONAL :: iunit
576 LOGICAL,
INTENT(IN),
OPTIONAL :: optbasis, allelectron, confine
578 INTEGER,
PARAMETER :: num_gto = 40
580 INTEGER :: i, ii, iw, k, l, ll, m, mb, mo, ngp, nn, &
581 nr, quadtype, relativistic, z
582 INTEGER,
DIMENSION(0:lmat) :: starti
583 INTEGER,
DIMENSION(0:lmat, 10) :: ncalc, ncore,
nelem
584 INTEGER,
DIMENSION(:),
POINTER :: econf
585 LOGICAL :: do_basopt, ecp_semi_local, monovalent
586 REAL(kind=
dp) :: al, aval, cc, cval, ear, rk, xx, zeff
587 REAL(kind=
dp),
DIMENSION(num_gto+2) :: results
602 NULLIFY (all_potential, gth_potential)
604 all_potential=all_potential, &
605 gth_potential=gth_potential, &
606 sgp_potential=sgp_potential, &
607 monovalent=monovalent)
609 IF (
PRESENT(iunit))
THEN
615 IF (
PRESENT(allelectron))
THEN
616 IF (allelectron)
THEN
617 NULLIFY (gth_potential)
623 IF (
PRESENT(optbasis))
THEN
627 cpassert(ngto <= num_gto)
629 IF (
ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential))
THEN
639 pp_calc=(
ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential)), &
641 relativistic=relativistic, &
645 ALLOCATE (potential, basis, integrals)
647 IF (
PRESENT(confine))
THEN
648 potential%confinement = confine
650 IF (
ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential))
THEN
651 potential%confinement = .true.
653 potential%confinement = .false.
657 potential%acon = 200._dp
658 potential%rcon = 4.0_dp
659 potential%scon = 8.0_dp
661 IF (
ASSOCIATED(gth_potential))
THEN
665 CALL set_atom(
atom, zcore=nint(zeff), potential=potential)
666 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
667 CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
668 IF (ecp_semi_local)
THEN
670 CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
671 potential%ecp_pot%symbol =
ptable(z)%symbol
674 CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
675 potential%sgp_pot%symbol =
ptable(z)%symbol
678 CALL set_atom(
atom, zcore=nint(zeff), potential=potential)
693 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
696 basis%eps_eig = 1.e-12_dp
699 basis%nprim = basis%nbas
700 m = maxval(basis%nbas)
701 ALLOCATE (basis%am(m, 0:
lmat))
704 DO i = 1, basis%nbas(l)
705 ll = i - 1 + starti(l)
706 basis%am(i, l) = aval*cval**(ll)
710 basis%geometrical = .true.
717 m = maxval(basis%nbas)
718 ALLOCATE (basis%bf(nr, m, 0:
lmat))
719 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
720 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
725 DO i = 1, basis%nbas(l)
728 rk = basis%grid%rad(k)
729 ear = exp(-al*basis%grid%rad(k)**2)
730 basis%bf(k, i, l) = rk**l*ear
731 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
732 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
733 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
741 atom%optimization%damping = 0.2_dp
742 atom%optimization%eps_scf = 1.e-6_dp
743 atom%optimization%eps_diis = 100._dp
744 atom%optimization%max_iter = 50
745 atom%optimization%n_diis = 5
753 ELSE IF (
ASSOCIATED(gth_potential))
THEN
756 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
760 DO l = 0, min(
lmat, ubound(
ptable(z)%e_conv, 1))
775 ncalc =
nelem - ncore
778 IF (qs_kind%ghost .OR. qs_kind%floating)
THEN
784 ALLOCATE (
atom%state)
786 atom%state%core = 0._dp
788 atom%state%occ = 0._dp
790 atom%state%occupation = 0._dp
791 atom%state%multiplicity = -1
795 IF (ncalc(l, i) > 0)
THEN
797 atom%state%occupation(l, k) = real(ncalc(l, i),
dp)
804 atom%state%maxl_calc =
atom%state%maxl_occ
805 atom%state%maxn_calc =
atom%state%maxn_occ
815 NULLIFY (integrals%tzora, integrals%hdkh)
820 mo = maxval(
atom%state%maxn_calc)
821 mb = maxval(
atom%basis%nbas)
832 density(i, 1) = xx*cc**i
833 density(i, 2) = results(2 + i)
837 density(1:ngto, 2) = results(1:ngto)
848 DEALLOCATE (potential, basis, integrals)
863 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: rtmat
865 INTEGER :: i, ii, ipgf, j, k, k1, k2, l, ll, m, n, &
866 ngp, nj, nn, nr, ns, nset, nsgf, &
867 quadtype, relativistic, z
868 INTEGER,
DIMENSION(0:lmat, 10) :: ncalc, ncore,
nelem
869 INTEGER,
DIMENSION(0:lmat, 100) :: set_index, shell_index
870 INTEGER,
DIMENSION(:),
POINTER :: lmax, lmin, npgf, nshell
871 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgf, last_sgf, ls
872 REAL(kind=
dp) :: al, alpha, ear, prefac, rk, zeff
873 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: omat
874 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: zet
875 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: gcc
884 IF (rel_control%rel_method ==
rel_none)
RETURN
886 NULLIFY (all_potential, orb_basis_set)
887 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, all_potential=all_potential)
889 cpassert(
ASSOCIATED(orb_basis_set))
891 IF (
ASSOCIATED(all_potential))
THEN
898 NULLIFY (
atom%xc_section)
899 NULLIFY (
atom%orbitals)
901 alpha = sqrt(all_potential%alpha_core_charge)
904 SELECT CASE (rel_control%rel_method)
906 cpabort(
"Unknown relativistic method for calculate_atomic_relkin")
908 SELECT CASE (rel_control%rel_DKH_order)
910 cpabort(
"Unknown DKH order for calculate_atomic_relkin")
921 SELECT CASE (rel_control%rel_zora_type)
923 cpabort(
"Unknown ZORA type for calculate_atomic_relkin")
925 cpabort(
"ZORA full not yet implemented in calculate_atomic_relkin")
936 relativistic=relativistic, &
940 ALLOCATE (potential, basis, integrals)
946 nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
947 first_sgf=first_sgf, last_sgf=last_sgf)
957 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
959 basis%eps_eig = 1.e-12_dp
967 DO j = lmin(i), min(lmax(i),
lmat)
968 basis%nprim(j) = basis%nprim(j) + npgf(i)
973 basis%nbas(l) = basis%nbas(l) + 1
977 shell_index(l, k) = j
982 nj = maxval(basis%nprim)
983 ns = maxval(basis%nbas)
984 ALLOCATE (basis%am(nj, 0:
lmat))
986 ALLOCATE (basis%cm(nj, ns, 0:
lmat))
992 IF (j >= lmin(i) .AND. j <= lmax(i))
THEN
994 basis%am(nj + ipgf, j) = zet(ipgf, i)
997 IF (ls(ii, i) == j)
THEN
1000 basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
1012 prefac = 2.0_dp*sqrt(
pi/
dfac(2*j + 1))
1013 DO ipgf = 1, basis%nprim(j)
1014 DO ii = 1, basis%nbas(j)
1015 basis%cm(ipgf, ii, j) = prefac*basis%cm(ipgf, ii, j)
1022 m = maxval(basis%nbas)
1023 ALLOCATE (basis%bf(nr, m, 0:
lmat))
1024 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
1025 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
1031 DO i = 1, basis%nprim(l)
1034 rk = basis%grid%rad(k)
1035 ear = exp(-al*basis%grid%rad(k)**2)
1036 DO j = 1, basis%nbas(l)
1037 basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
1038 basis%dbf(k, j, l) = basis%dbf(k, j, l) &
1039 + (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
1040 basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
1041 (real(l*(l - 1),
dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1,
dp)* &
1042 rk**(l) + 4._dp*al*rk**(l + 2))*ear*basis%cm(i, j, l)
1051 atom%optimization%damping = 0.2_dp
1052 atom%optimization%eps_scf = 1.e-6_dp
1053 atom%optimization%eps_diis = 100._dp
1054 atom%optimization%max_iter = 50
1055 atom%optimization%n_diis = 5
1061 DO l = 0, min(
lmat, ubound(
ptable(z)%e_conv, 1))
1076 ncalc =
nelem - ncore
1078 IF (qs_kind%ghost .OR. qs_kind%floating)
THEN
1084 ALLOCATE (
atom%state)
1086 atom%state%core = 0._dp
1088 atom%state%occ = 0._dp
1090 atom%state%occupation = 0._dp
1091 atom%state%multiplicity = -1
1095 IF (ncalc(l, i) > 0)
THEN
1097 atom%state%occupation(l, k) = real(ncalc(l, i),
dp)
1104 atom%state%maxl_calc =
atom%state%maxl_occ
1105 atom%state%maxn_calc =
atom%state%maxn_occ
1113 NULLIFY (integrals%tzora, integrals%hdkh)
1119 integrals%core = 0.0_dp
1123 ALLOCATE (omat(m, m))
1125 CALL sg_erfc(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
1126 integrals%core(1:n, 1:n, l) = matmul(transpose(basis%cm(1:m, 1:n, l)), &
1127 matmul(omat(1:m, 1:m), basis%cm(1:m, 1:n, l)))
1133 IF (
ASSOCIATED(rtmat))
THEN
1136 ALLOCATE (rtmat(nsgf, nsgf))
1140 DO k1 = 1, basis%nbas(l)
1141 DO k2 = 1, basis%nbas(l)
1142 i = first_sgf(shell_index(l, k1), set_index(l, k1))
1143 j = first_sgf(shell_index(l, k2), set_index(l, k2))
1144 SELECT CASE (
atom%relativistic)
1146 cpabort(
"Unknown relativistic type for calculate_atomic_relkin")
1149 rtmat(i + m, j + m) = integrals%tzora(k1, k2, l)
1153 rtmat(i + m, j + m) = integrals%hdkh(k1, k2, l) - integrals%kin(k1, k2, l) + &
1154 atom%zcore*integrals%core(k1, k2, l)
1162 rtmat(k1, k2) = 0.5_dp*(rtmat(k1, k2) + rtmat(k2, k1))
1163 rtmat(k2, k1) = rtmat(k1, k2)
1175 DEALLOCATE (potential, basis, integrals)
1179 IF (
ASSOCIATED(rtmat))
THEN
1197 INTEGER :: i, j, l, lm, n, ne, nexp_lpot, nexp_lsd, &
1199 INTEGER,
DIMENSION(:),
POINTER :: nct_lpot, nct_lsd, nct_nlcc, nppnl, &
1201 LOGICAL :: lpot_present, lsd_present, nlcc_present, &
1203 REAL(kind=
dp) :: ac, zeff
1204 REAL(kind=
dp),
DIMENSION(:),
POINTER :: alpha_lpot, alpha_lsd, alpha_nlcc, ap, ce
1205 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: cval_lpot, cval_lsd, cval_nlcc
1206 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: hp, kp
1210 elec_conf=ppeconf, &
1211 alpha_core_charge=ac, &
1220 gth_atompot%zion = zeff
1221 gth_atompot%rc = sqrt(0.5_dp/ac)
1222 gth_atompot%ncl = ne
1223 gth_atompot%cl(:) = 0._dp
1224 IF (ac > 0._dp)
THEN
1226 gth_atompot%cl(i) = ce(i)/(2._dp*ac)**(i - 1)
1230 gth_atompot%lpotextended = .false.
1231 gth_atompot%lsdpot = .false.
1232 gth_atompot%nlcc = .false.
1233 gth_atompot%nexp_lpot = 0
1234 gth_atompot%nexp_lsd = 0
1235 gth_atompot%nexp_nlcc = 0
1237 lpot_present=lpot_present, &
1238 lsd_present=lsd_present, &
1239 nlcc_present=nlcc_present)
1240 IF (lpot_present)
THEN
1242 nexp_lpot=nexp_lpot, &
1243 alpha_lpot=alpha_lpot, &
1244 nct_lpot=nct_lpot, &
1245 cval_lpot=cval_lpot)
1246 gth_atompot%lpotextended = .true.
1247 gth_atompot%nexp_lpot = nexp_lpot
1248 gth_atompot%alpha_lpot(1:nexp_lpot) = sqrt(0.5_dp/alpha_lpot(1:nexp_lpot))
1249 gth_atompot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot)
1253 gth_atompot%cval_lpot(i, j) = cval_lpot(i, j)/(2._dp*ac)**(i - 1)
1257 IF (lsd_present)
THEN
1259 nexp_lsd=nexp_lsd, &
1260 alpha_lsd=alpha_lsd, &
1263 gth_atompot%lsdpot = .true.
1264 gth_atompot%nexp_lsd = nexp_lsd
1265 gth_atompot%alpha_lsd(1:nexp_lsd) = sqrt(0.5_dp/alpha_lsd(1:nexp_lsd))
1266 gth_atompot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd)
1270 gth_atompot%cval_lsd(i, j) = cval_lsd(i, j)/(2._dp*ac)**(i - 1)
1276 gth_atompot%nl(:) = 0
1277 gth_atompot%rcnl(:) = 0._dp
1278 gth_atompot%hnl(:, :, :) = 0._dp
1281 gth_atompot%nl(l) = n
1282 gth_atompot%rcnl(l) = sqrt(0.5_dp/ap(l))
1283 gth_atompot%hnl(1:n, 1:n, l) = hp(1:n, 1:n, l)
1288 gth_atompot%soc = soc_present
1289 gth_atompot%knl = 0.0_dp
1290 IF (soc_present)
THEN
1293 gth_atompot%knl(1:n, 1:n, l) = kp(1:n, 1:n, l)
1297 IF (nlcc_present)
THEN
1299 nexp_nlcc=nexp_nlcc, &
1300 alpha_nlcc=alpha_nlcc, &
1301 nct_nlcc=nct_nlcc, &
1302 cval_nlcc=cval_nlcc)
1303 gth_atompot%nlcc = .true.
1304 gth_atompot%nexp_nlcc = nexp_nlcc
1305 gth_atompot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc)
1306 gth_atompot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc)
1307 gth_atompot%cval_nlcc(1:4, 1:nexp_nlcc) = cval_nlcc(1:4, 1:nexp_nlcc)
1317 SUBROUTINE sgp_potential_conversion(sgp_potential, sgp_atompot)
1322 INTEGER,
DIMENSION(:),
POINTER :: ppeconf
1323 LOGICAL :: nlcc_present
1324 REAL(kind=
dp) :: ac, zeff
1325 REAL(kind=
dp),
DIMENSION(:),
POINTER :: ap, ce
1326 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: hhp
1327 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: ccp
1330 name=sgp_atompot%pname, &
1332 elec_conf=ppeconf, &
1333 alpha_core_charge=ac)
1334 sgp_atompot%zion = zeff
1335 sgp_atompot%ac_local = ac
1336 sgp_atompot%econf(0:3) = ppeconf(0:3)
1338 is_nonlocal=sgp_atompot%is_nonlocal, &
1339 n_nonlocal=n, a_nonlocal=ap, h_nonlocal=hhp, c_nonlocal=ccp)
1341 sgp_atompot%has_nonlocal = any(sgp_atompot%is_nonlocal)
1342 sgp_atompot%lmax = lm
1343 IF (sgp_atompot%has_nonlocal)
THEN
1344 cpassert(n <=
SIZE(sgp_atompot%a_nonlocal))
1345 sgp_atompot%n_nonlocal = n
1346 sgp_atompot%a_nonlocal(1:n) = ap(1:n)
1347 sgp_atompot%h_nonlocal(1:n, 0:lm) = hhp(1:n, 0:lm)
1348 sgp_atompot%c_nonlocal(1:n, 1:n, 0:lm) = ccp(1:n, 1:n, 0:lm)
1351 CALL get_potential(sgp_potential, n_local=n, a_local=ap, c_local=ce)
1352 cpassert(n <=
SIZE(sgp_atompot%a_local))
1353 sgp_atompot%n_local = n
1354 sgp_atompot%a_local(1:n) = ap(1:n)
1355 sgp_atompot%c_local(1:n) = ce(1:n)
1358 n_nlcc=n, a_nlcc=ap, c_nlcc=ce)
1359 IF (nlcc_present)
THEN
1360 sgp_atompot%has_nlcc = .true.
1361 cpassert(n <=
SIZE(sgp_atompot%a_nlcc))
1362 sgp_atompot%n_nlcc = n
1363 sgp_atompot%a_nlcc(1:n) = ap(1:n)
1364 sgp_atompot%c_nlcc(1:n) = ce(1:n)
1366 sgp_atompot%has_nlcc = .false.
1369 END SUBROUTINE sgp_potential_conversion
1376 SUBROUTINE ecp_potential_conversion(sgp_potential, ecp_atompot)
1380 INTEGER,
DIMENSION(:),
POINTER :: ppeconf
1381 LOGICAL :: ecp_local, ecp_semi_local
1382 REAL(kind=
dp) :: zeff
1384 CALL get_potential(sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
1385 cpassert(ecp_semi_local .AND. ecp_local)
1387 name=ecp_atompot%pname, &
1390 ecp_atompot%zion = zeff
1391 ecp_atompot%econf(0:3) = ppeconf(0:3)
1394 CALL get_potential(sgp_potential, nloc=ecp_atompot%nloc, nrloc=ecp_atompot%nrloc, &
1395 aloc=ecp_atompot%aloc, bloc=ecp_atompot%bloc)
1397 CALL get_potential(sgp_potential, npot=ecp_atompot%npot, nrpot=ecp_atompot%nrpot, &
1398 apot=ecp_atompot%apot, bpot=ecp_atompot%bpot)
1400 END SUBROUTINE ecp_potential_conversion
subroutine, public sg_erfc(umat, l, a, pa, pb)
...
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
routines that fit parameters for /from atomic calculations
subroutine, public atom_fit_density(atom, num_gto, norder, iunit, agto, powell_section, results)
Fit the atomic electron density using a geometrical Gaussian basis set.
calculate the orbitals for a given atomic kind type
subroutine, public calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
...
subroutine, public calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, optbasis, allelectron, confine)
...
subroutine, public gth_potential_conversion(gth_potential, gth_atompot)
...
subroutine, public calculate_atomic_orbitals(atomic_kind, qs_kind, agrid, iunit, pmat, fmat, density, wavefunction, wfninfo, confine, xc_section, nocc, which_l, which_n, proj_shell_charge, ao_coef)
...
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_ppint_setup(integrals, basis, potential)
...
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
subroutine, public set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid, cp2k_norm)
...
Define the atom type and its sub types.
subroutine, public create_atom_type(atom)
...
integer, parameter, public cgto_basis
integer, parameter, public gto_basis
subroutine, public release_atom_type(atom)
...
subroutine, public release_atom_potential(potential)
...
integer, parameter, public lmat
subroutine, public set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
...
subroutine, public release_atom_basis(basis)
...
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
subroutine, public clementi_geobas(zval, cval, aval, ngto, ival)
...
Some basic routines for atomic calculations.
pure integer function, dimension(0:lmat), public get_maxn_occ(occupation)
Return the maximum principal quantum number of occupied orbitals.
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
pure integer function, public get_maxl_occ(occupation)
Return the maximum orbital quantum number of occupied orbitals.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition of the atomic potential types.
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
integer, parameter, public nelem
Definition of physical constants:
real(kind=dp), parameter, public bohr
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Define the quickstep kind type and their sub types.
subroutine, public set_pseudo_state(econf, z, ncalc, ncore, nelem)
...
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public init_atom_electronic_state(atomic_kind, qs_kind, ncalc, ncore, nelem, edelta, rks)
...
parameters that control a relativistic calculation
Provides all information about a basis set.
Provides all information about a pseudopotential.
Holds atomic orbitals and energies.
Provides all information about an atomic kind.
Provides all information about an atomic kind.
Provides all information about a quickstep kind.
contains the parameters needed by a relativistic calculation