59#include "./base/base_uses.f90"
65 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'atom_types'
68 INTEGER,
PARAMETER ::
lmat = 5
75 INTEGER,
PARAMETER :: nmax = 25
81 INTEGER,
DIMENSION(0:lmat) :: nbas = 0
82 INTEGER,
DIMENSION(0:lmat) :: nprim = 0
83 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: am => null()
84 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: cm => null()
85 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: as => null()
86 INTEGER,
DIMENSION(:, :),
POINTER :: ns => null()
87 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: bf => null()
88 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: dbf => null()
89 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: ddbf => null()
90 REAL(kind=
dp) :: eps_eig = 0.0_dp
92 LOGICAL :: geometrical = .false.
93 REAL(kind=
dp) :: aval = 0.0_dp, cval = 0.0_dp
94 INTEGER,
DIMENSION(0:lmat) :: start = 0
100 CHARACTER(LEN=2) :: symbol =
""
101 CHARACTER(LEN=default_string_length) :: pname =
""
102 INTEGER,
DIMENSION(0:lmat) :: econf = 0
103 REAL(
dp) :: zion = 0.0_dp
104 REAL(
dp) :: rc = 0.0_dp
106 REAL(
dp),
DIMENSION(5) :: cl = 0.0_dp
107 INTEGER,
DIMENSION(0:lmat) :: nl = 0
108 REAL(
dp),
DIMENSION(0:lmat) :: rcnl = 0.0_dp
109 REAL(
dp),
DIMENSION(4, 4, 0:lmat) :: hnl = 0.0_dp
111 LOGICAL :: soc = .false.
112 REAL(
dp),
DIMENSION(4, 4, 0:lmat) :: knl = 0.0_dp
115 LOGICAL :: nlcc = .false.
116 INTEGER :: nexp_nlcc = 0
117 REAL(kind=
dp),
DIMENSION(10) :: alpha_nlcc = 0.0_dp
118 INTEGER,
DIMENSION(10) :: nct_nlcc = 0
119 REAL(kind=
dp),
DIMENSION(4, 10) :: cval_nlcc = 0.0_dp
121 LOGICAL :: lsdpot = .false.
122 INTEGER :: nexp_lsd = 0
123 REAL(kind=
dp),
DIMENSION(10) :: alpha_lsd = 0.0_dp
124 INTEGER,
DIMENSION(10) :: nct_lsd = 0
125 REAL(kind=
dp),
DIMENSION(4, 10) :: cval_lsd = 0.0_dp
127 LOGICAL :: lpotextended = .false.
128 INTEGER :: nexp_lpot = 0
129 REAL(kind=
dp),
DIMENSION(10) :: alpha_lpot = 0.0_dp
130 INTEGER,
DIMENSION(10) :: nct_lpot = 0
131 REAL(kind=
dp),
DIMENSION(4, 10) :: cval_lpot = 0.0_dp
135 CHARACTER(LEN=2) :: symbol =
""
136 CHARACTER(LEN=default_string_length) :: pname =
""
137 INTEGER,
DIMENSION(0:lmat) :: econf = 0
138 REAL(
dp) :: zion = 0.0_dp
141 INTEGER,
DIMENSION(1:15) :: nrloc = 0
142 REAL(
dp),
DIMENSION(1:15) :: aloc = 0.0_dp
143 REAL(
dp),
DIMENSION(1:15) :: bloc = 0.0_dp
144 INTEGER,
DIMENSION(0:10) :: npot = 0
145 INTEGER,
DIMENSION(1:15, 0:10) :: nrpot = 0
146 REAL(
dp),
DIMENSION(1:15, 0:10) :: apot = 0.0_dp
147 REAL(
dp),
DIMENSION(1:15, 0:10) :: bpot = 0.0_dp
151 CHARACTER(LEN=2) :: symbol =
""
152 CHARACTER(LEN=default_string_length) :: pname =
""
153 INTEGER,
DIMENSION(0:lmat) :: econf = 0
154 REAL(
dp) :: zion = 0.0_dp
156 LOGICAL :: has_nonlocal = .false.
157 INTEGER :: n_nonlocal = 0
158 LOGICAL,
DIMENSION(0:5) :: is_nonlocal = .false.
159 REAL(kind=
dp),
DIMENSION(nmax) :: a_nonlocal = 0.0_dp
160 REAL(kind=
dp),
DIMENSION(nmax, 0:lmat) :: h_nonlocal = 0.0_dp
161 REAL(kind=
dp),
DIMENSION(nmax, nmax, 0:lmat) :: c_nonlocal = 0.0_dp
162 INTEGER :: n_local = 0
163 REAL(kind=
dp) :: ac_local = 0.0_dp
164 REAL(kind=
dp),
DIMENSION(nmax) :: a_local = 0.0_dp
165 REAL(kind=
dp),
DIMENSION(nmax) :: c_local = 0.0_dp
166 LOGICAL :: has_nlcc = .false.
167 INTEGER :: n_nlcc = 0
168 REAL(kind=
dp),
DIMENSION(nmax) :: a_nlcc = 0.0_dp
169 REAL(kind=
dp),
DIMENSION(nmax) :: c_nlcc = 0.0_dp
173 INTEGER :: ppot_type = 0
174 LOGICAL :: confinement = .false.
175 INTEGER :: conf_type = 0
176 REAL(
dp) :: acon = 0.0_dp
177 REAL(
dp) :: rcon = 0.0_dp
178 REAL(
dp) :: scon = 0.0_dp
181 TYPE(atom_upfpot_type) :: upf_pot = atom_upfpot_type()
189 REAL(kind=
dp) :: scale_coulomb = 0.0_dp
190 REAL(kind=
dp) :: scale_longrange = 0.0_dp
191 REAL(kind=
dp) :: omega = 0.0_dp
192 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: kernel
193 LOGICAL :: do_gh = .false.
200 REAL(kind=
dp),
DIMENSION(0:lmat, 10) :: occ = 0.0_dp
201 REAL(kind=
dp),
DIMENSION(0:lmat, 10) :: core = 0.0_dp
202 REAL(kind=
dp),
DIMENSION(0:lmat, 10) :: occupation = 0.0_dp
203 INTEGER :: maxl_occ = 0
204 INTEGER,
DIMENSION(0:lmat) :: maxn_occ = 0
205 INTEGER :: maxl_calc = 0
206 INTEGER,
DIMENSION(0:lmat) :: maxn_calc = 0
207 INTEGER :: multiplicity = 0
208 REAL(kind=
dp),
DIMENSION(0:lmat, 10) :: occa = 0.0_dp, occb = 0.0_dp
214 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: int => null()
218 INTEGER :: status = 0
219 INTEGER :: ppstat = 0
220 LOGICAL :: eri_coulomb = .false.
221 LOGICAL :: eri_exchange = .false.
222 LOGICAL :: all_nu = .false.
223 INTEGER,
DIMENSION(0:lmat) :: n = 0, nne = 0
224 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: ovlp => null(), kin => null(), core => null(), clsd => null()
225 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: utrans => null(), uptrans => null()
226 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: hnl => null()
227 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: conf => null()
228 TYPE(
eri),
DIMENSION(100) :: ceri =
eri()
229 TYPE(
eri),
DIMENSION(100) :: eeri =
eri()
230 INTEGER :: dkhstat = 0
231 INTEGER :: zorastat = 0
232 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: tzora => null()
233 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: hdkh => null()
239 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: wfn => null(), wfna => null(), wfnb => null()
240 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: pmat => null(), pmata => null(), pmatb => null()
241 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ener => null(), enera => null(), enerb => null()
242 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: refene => null(), refchg => null(), refnod => null()
243 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: wrefene => null(), wrefchg => null(), wrefnod => null()
244 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: crefene => null(), crefchg => null(), crefnod => null()
245 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: wpsir0 => null(), tpsir0 => null()
246 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: rcmax => null()
247 CHARACTER(LEN=2),
DIMENSION(:, :, :),
POINTER :: reftype => null()
253 INTEGER,
DIMENSION(0:lmat) :: n = 0
254 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: op => null()
260 REAL(kind=
dp),
DIMENSION(:),
POINTER :: op => null()
266 TYPE atom_energy_type
267 REAL(kind=
dp) :: etot = 0.0_dp
268 REAL(kind=
dp) :: eband = 0.0_dp
269 REAL(kind=
dp) :: ekin = 0.0_dp
270 REAL(kind=
dp) :: epot = 0.0_dp
271 REAL(kind=
dp) :: ecore = 0.0_dp
272 REAL(kind=
dp) :: elsd = 0.0_dp
273 REAL(kind=
dp) :: epseudo = 0.0_dp
274 REAL(kind=
dp) :: eploc = 0.0_dp
275 REAL(kind=
dp) :: epnl = 0.0_dp
276 REAL(kind=
dp) :: exc = 0.0_dp
277 REAL(kind=
dp) :: ecoulomb = 0.0_dp
278 REAL(kind=
dp) :: eexchange = 0.0_dp
279 REAL(kind=
dp) :: econfinement = 0.0_dp
280 END TYPE atom_energy_type
285 REAL(kind=
dp) :: damping = 0.0_dp
286 REAL(kind=
dp) :: eps_scf = 0.0_dp
287 REAL(kind=
dp) :: eps_diis = 0.0_dp
288 INTEGER :: max_iter = 0
289 INTEGER :: n_diis = 0
297 LOGICAL :: pp_calc = .false.
299 LOGICAL :: do_zmp = .false., doread = .false., read_vxc = .false., dm = .false.
300 CHARACTER(LEN=default_string_length) :: ext_file =
"", ext_vxc_file =
"", &
301 zmp_restart_file =
""
308 REAL(kind=
dp) :: lambda = 0.0_dp
309 REAL(kind=
dp) :: rho_diff_integral = 0.0_dp
310 REAL(kind=
dp) :: weight = 0.0_dp, zmpgrid_tol = 0.0_dp, zmpvxcgrid_tol = 0.0_dp
317 TYPE(atom_energy_type) :: energy = atom_energy_type()
347 MODULE PROCEDURE read_ecp_potential_file, &
348 read_ecp_potential_files
386 INTEGER,
INTENT(IN) :: zval
387 CHARACTER(LEN=2) :: btyp
389 CHARACTER(LEN=*),
PARAMETER :: routinen =
'init_atom_basis'
390 INTEGER,
PARAMETER :: nua = 40, nup = 16
391 REAL(kind=
dp),
DIMENSION(nua),
PARAMETER :: ugbs = [0.007299_dp, 0.013705_dp, 0.025733_dp, &
392 0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
393 3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
394 174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
395 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
396 94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
397 2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
398 27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
401 CHARACTER(LEN=default_path_length) :: basis_fn
402 CHARACTER(LEN=default_string_length) :: basis_name
403 INTEGER :: basistype, handle, i, j, k, l, ll, m, &
404 ngp, nl, nr, nu, quadtype
405 INTEGER,
DIMENSION(0:lmat) :: starti
406 INTEGER,
DIMENSION(:),
POINTER :: nqm, num_gto, num_slater, sindex
407 REAL(kind=
dp) :: al, amax, aval, cval, ear, pf, rk
408 REAL(kind=
dp),
DIMENSION(:),
POINTER :: expo
411 CALL timeset(routinen, handle)
418 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
425 cpabort(
"The number of radial grid points must be greater than zero.")
429 basis%geometrical = .false.
436 SELECT CASE (basistype)
441 IF (num_gto(1) < 1)
THEN
443 IF (btyp ==
"AE")
THEN
445 ELSE IF (btyp ==
"PP")
THEN
452 ALLOCATE (basis%am(nu, 0:
lmat))
454 basis%am(1:nu, i) = ugbs(1:nu)
458 DO i = 1,
SIZE(num_gto)
459 basis%nbas(i - 1) = num_gto(i)
461 basis%nprim = basis%nbas
462 m = maxval(basis%nbas)
463 ALLOCATE (basis%am(m, 0:
lmat))
466 IF (basis%nbas(l) > 0)
THEN
478 cpabort(
"Invalid angular quantum number l found for Gaussian basis set")
480 cpassert(
SIZE(expo) >= basis%nbas(l))
481 DO i = 1, basis%nbas(l)
482 basis%am(i, l) = expo(i)
489 m = maxval(basis%nbas)
490 ALLOCATE (basis%bf(nr, m, 0:
lmat))
491 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
492 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
497 DO i = 1, basis%nbas(l)
500 rk = basis%grid%rad(k)
501 ear = exp(-al*basis%grid%rad(k)**2)
502 basis%bf(k, i, l) = rk**l*ear
503 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
504 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
505 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
513 IF (num_gto(1) < 1)
THEN
514 IF (btyp ==
"AE")
THEN
517 ELSE IF (btyp ==
"PP")
THEN
520 ELSE IF (btyp ==
"AA")
THEN
522 amax = cval**(basis%nbas(0) - 1)
523 basis%nbas(0) = nint((log(amax)/log(1.6_dp)))
526 basis%nbas(1) = basis%nbas(0) - 4
527 basis%nbas(2) = basis%nbas(0) - 8
528 basis%nbas(3) = basis%nbas(0) - 12
529 IF (
lmat > 3) basis%nbas(4:
lmat) = 0
530 ELSE IF (btyp ==
"AP")
THEN
533 basis%nbas = nint((log(amax)/log(1.6_dp)))
540 basis%nprim = basis%nbas
543 DO i = 1,
SIZE(num_gto)
544 basis%nbas(i - 1) = num_gto(i)
546 basis%nprim = basis%nbas
550 DO i = 1,
SIZE(sindex)
551 starti(i - 1) = sindex(i)
552 cpassert(sindex(i) >= 0)
557 m = maxval(basis%nbas)
558 ALLOCATE (basis%am(m, 0:
lmat))
561 DO i = 1, basis%nbas(l)
562 ll = i - 1 + starti(l)
563 basis%am(i, l) = aval*cval**(ll)
567 basis%geometrical = .true.
574 m = maxval(basis%nbas)
575 ALLOCATE (basis%bf(nr, m, 0:
lmat))
576 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
577 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
582 DO i = 1, basis%nbas(l)
585 rk = basis%grid%rad(k)
586 ear = exp(-al*basis%grid%rad(k)**2)
587 basis%bf(k, i, l) = rk**l*ear
588 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
589 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
590 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
599 CALL read_basis_set(
ptable(zval)%symbol, basis, basis_name, basis_fn, &
604 m = maxval(basis%nbas)
605 ALLOCATE (basis%bf(nr, m, 0:
lmat))
606 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
607 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
612 DO i = 1, basis%nprim(l)
615 rk = basis%grid%rad(k)
616 ear = exp(-al*basis%grid%rad(k)**2)
617 DO j = 1, basis%nbas(l)
618 basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
619 basis%dbf(k, j, l) = basis%dbf(k, j, l) &
620 + (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
621 basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
622 (real(l*(l - 1),
dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))* &
623 ear*basis%cm(i, j, l)
632 IF (num_slater(1) < 1)
THEN
633 cpabort(
"Invalid number (less than 1) Slater-type functions found.")
636 DO i = 1,
SIZE(num_slater)
637 basis%nbas(i - 1) = num_slater(i)
639 basis%nprim = basis%nbas
640 m = maxval(basis%nbas)
641 ALLOCATE (basis%as(m, 0:
lmat), basis%ns(m, 0:
lmat))
645 IF (basis%nbas(l) > 0)
THEN
657 cpabort(
"Invalid angular quantum number l found for Slater basis set")
659 cpassert(
SIZE(expo) >= basis%nbas(l))
660 DO i = 1, basis%nbas(l)
661 basis%as(i, l) = expo(i)
674 cpabort(
"Invalid angular quantum number l found for Slater basis set")
676 cpassert(
SIZE(nqm) >= basis%nbas(l))
677 DO i = 1, basis%nbas(l)
678 basis%ns(i, l) = nqm(i)
685 m = maxval(basis%nbas)
686 ALLOCATE (basis%bf(nr, m, 0:
lmat))
687 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
688 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
693 DO i = 1, basis%nbas(l)
696 pf = (2._dp*al)**nl*sqrt(2._dp*al/
fac(2*nl))
698 rk = basis%grid%rad(k)
699 ear = rk**(nl - 1)*exp(-al*rk)
700 basis%bf(k, i, l) = pf*ear
701 basis%dbf(k, i, l) = pf*(real(nl - 1,
dp)/rk - al)*ear
702 basis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1),
dp)/rk/rk &
703 - al*real(2*(nl - 1),
dp)/rk + al*al)*ear
709 cpabort(
"Numerical basis set type not yet implemented.")
711 cpabort(
"Unknown basis set type specified. Check the code!")
714 CALL timestop(handle)
725 CHARACTER(LEN=*),
PARAMETER :: routinen =
'init_atom_basis_default_pp'
726 INTEGER,
PARAMETER :: nua = 40, nup = 20
727 REAL(kind=
dp),
DIMENSION(nua),
PARAMETER :: ugbs = [0.007299_dp, 0.013705_dp, 0.025733_dp, &
728 0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
729 3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
730 174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
731 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
732 94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
733 2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
734 27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
737 INTEGER :: handle, i, k, l, m, ngp, nr, nu, quadtype
738 REAL(kind=
dp) :: al, ear, rk
740 CALL timeset(routinen, handle)
742 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
751 basis%geometrical = .false.
755 basis%eps_eig = 1.e-12_dp
761 ALLOCATE (basis%am(nu, 0:
lmat))
763 basis%am(1:nu, i) = ugbs(1:nu)
767 m = maxval(basis%nbas)
768 ALLOCATE (basis%bf(nr, m, 0:
lmat))
769 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
770 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
775 DO i = 1, basis%nbas(l)
778 rk = basis%grid%rad(k)
779 ear = exp(-al*basis%grid%rad(k)**2)
780 basis%bf(k, i, l) = rk**l*ear
781 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
782 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
783 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
788 CALL timestop(handle)
802 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: r, rab
804 INTEGER :: i, j, k, l, m, n1, n2, n3, ngp, nl, nr
805 REAL(kind=
dp) :: al, ear, pf, rk
807 NULLIFY (gbasis%am, gbasis%cm, gbasis%as, gbasis%ns, gbasis%bf, gbasis%dbf, gbasis%ddbf)
810 gbasis%basis_type = basis%basis_type
811 gbasis%nbas(0:
lmat) = basis%nbas(0:
lmat)
812 gbasis%nprim(0:
lmat) = basis%nprim(0:
lmat)
813 IF (
ASSOCIATED(basis%am))
THEN
814 n1 =
SIZE(basis%am, 1)
815 n2 =
SIZE(basis%am, 2)
816 ALLOCATE (gbasis%am(n1, 0:n2 - 1))
819 IF (
ASSOCIATED(basis%cm))
THEN
820 n1 =
SIZE(basis%cm, 1)
821 n2 =
SIZE(basis%cm, 2)
822 n3 =
SIZE(basis%cm, 3)
823 ALLOCATE (gbasis%cm(n1, n2, 0:n3 - 1))
826 IF (
ASSOCIATED(basis%as))
THEN
827 n1 =
SIZE(basis%as, 1)
828 n2 =
SIZE(basis%as, 2)
829 ALLOCATE (gbasis%as(n1, 0:n2 - 1))
832 IF (
ASSOCIATED(basis%ns))
THEN
833 n1 =
SIZE(basis%ns, 1)
834 n2 =
SIZE(basis%ns, 2)
835 ALLOCATE (gbasis%ns(n1, 0:n2 - 1))
838 gbasis%eps_eig = basis%eps_eig
839 gbasis%geometrical = basis%geometrical
840 gbasis%aval = basis%aval
841 gbasis%cval = basis%cval
842 gbasis%start(0:
lmat) = basis%start(0:
lmat)
846 NULLIFY (gbasis%grid)
850 cpabort(
"The number of radial grid points must be greater than zero.")
852 cpassert(
SIZE(rab) == ngp)
854 ALLOCATE (gbasis%grid%rad(ngp), gbasis%grid%rad2(ngp), gbasis%grid%wr(ngp))
856 gbasis%grid%rad(:) = r(:)
857 gbasis%grid%rad2(:) = r(:)*r(:)
858 gbasis%grid%wr(:) = rab(:)*gbasis%grid%rad2(:)
862 m = maxval(gbasis%nbas)
863 ALLOCATE (gbasis%bf(nr, m, 0:
lmat))
864 ALLOCATE (gbasis%dbf(nr, m, 0:
lmat))
865 ALLOCATE (gbasis%ddbf(nr, m, 0:
lmat))
870 SELECT CASE (gbasis%basis_type)
874 DO i = 1, gbasis%nbas(l)
877 rk = gbasis%grid%rad(k)
878 ear = exp(-al*gbasis%grid%rad(k)**2)
879 gbasis%bf(k, i, l) = rk**l*ear
880 gbasis%dbf(k, i, l) = (real(l,
dp)*rk**max(0, l - 1) - 2._dp*al*rk**(l + 1))*ear
881 gbasis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**max(0, l - 2) - &
882 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
888 DO i = 1, gbasis%nprim(l)
891 rk = gbasis%grid%rad(k)
892 ear = exp(-al*gbasis%grid%rad(k)**2)
893 DO j = 1, gbasis%nbas(l)
894 gbasis%bf(k, j, l) = gbasis%bf(k, j, l) + rk**l*ear*gbasis%cm(i, j, l)
895 gbasis%dbf(k, j, l) = gbasis%dbf(k, j, l) &
896 + (real(l,
dp)*rk**max(0, l - 1) - 2._dp*al*rk**(l + 1))*ear*gbasis%cm(i, j, l)
897 gbasis%ddbf(k, j, l) = gbasis%ddbf(k, j, l) + &
898 (real(l*(l - 1),
dp)*rk**max(0, l - 2) - 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))* &
899 ear*gbasis%cm(i, j, l)
906 DO i = 1, gbasis%nbas(l)
909 pf = (2._dp*al)**nl*sqrt(2._dp*al/
fac(2*nl))
911 rk = gbasis%grid%rad(k)
912 ear = rk**(nl - 1)*exp(-al*rk)
913 gbasis%bf(k, i, l) = pf*ear
914 gbasis%dbf(k, i, l) = pf*(real(nl - 1,
dp)/rk - al)*ear
915 gbasis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1),
dp)/rk/rk &
916 - al*real(2*(nl - 1),
dp)/rk + al*al)*ear
922 cpabort(
"Numerical basis set type not yet implemented.")
924 cpabort(
"Unknown basis set type specified. Check the code!")
936 IF (
ASSOCIATED(basis%am))
THEN
937 DEALLOCATE (basis%am)
939 IF (
ASSOCIATED(basis%cm))
THEN
940 DEALLOCATE (basis%cm)
942 IF (
ASSOCIATED(basis%as))
THEN
943 DEALLOCATE (basis%as)
945 IF (
ASSOCIATED(basis%ns))
THEN
946 DEALLOCATE (basis%ns)
948 IF (
ASSOCIATED(basis%bf))
THEN
949 DEALLOCATE (basis%bf)
951 IF (
ASSOCIATED(basis%dbf))
THEN
952 DEALLOCATE (basis%dbf)
954 IF (
ASSOCIATED(basis%ddbf))
THEN
955 DEALLOCATE (basis%ddbf)
970 cpassert(.NOT.
ASSOCIATED(
atom))
974 NULLIFY (
atom%zmp_section)
975 NULLIFY (
atom%xc_section)
977 atom%do_zmp = .false.
978 atom%doread = .false.
979 atom%read_vxc = .false.
981 atom%hfx_pot%scale_coulomb = 0.0_dp
982 atom%hfx_pot%scale_longrange = 0.0_dp
983 atom%hfx_pot%omega = 0.0_dp
994 cpassert(
ASSOCIATED(
atom))
997 NULLIFY (
atom%integrals)
998 IF (
ASSOCIATED(
atom%state))
THEN
999 DEALLOCATE (
atom%state)
1001 IF (
ASSOCIATED(
atom%orbitals))
THEN
1031 SUBROUTINE set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, &
1032 read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
1039 INTEGER,
INTENT(IN),
OPTIONAL :: zcore
1040 LOGICAL,
INTENT(IN),
OPTIONAL :: pp_calc, do_zmp, doread, read_vxc
1041 INTEGER,
INTENT(IN),
OPTIONAL :: method_type, relativistic, &
1042 coulomb_integral_type, &
1043 exchange_integral_type
1046 cpassert(
ASSOCIATED(
atom))
1048 IF (
PRESENT(basis))
atom%basis => basis
1049 IF (
PRESENT(state))
atom%state => state
1050 IF (
PRESENT(integrals))
atom%integrals => integrals
1051 IF (
PRESENT(orbitals))
atom%orbitals => orbitals
1052 IF (
PRESENT(potential))
atom%potential => potential
1053 IF (
PRESENT(zcore))
atom%zcore = zcore
1054 IF (
PRESENT(pp_calc))
atom%pp_calc = pp_calc
1056 IF (
PRESENT(do_zmp))
atom%do_zmp = do_zmp
1057 IF (
PRESENT(doread))
atom%doread = doread
1058 IF (
PRESENT(read_vxc))
atom%read_vxc = read_vxc
1060 IF (
PRESENT(method_type))
atom%method_type = method_type
1061 IF (
PRESENT(relativistic))
atom%relativistic = relativistic
1062 IF (
PRESENT(coulomb_integral_type))
atom%coulomb_integral_type = coulomb_integral_type
1063 IF (
PRESENT(exchange_integral_type))
atom%exchange_integral_type = exchange_integral_type
1065 IF (
PRESENT(fmat))
THEN
1080 INTEGER,
INTENT(IN) :: mbas, mo
1082 cpassert(.NOT.
ASSOCIATED(orbs))
1086 ALLOCATE (orbs%wfn(mbas, mo, 0:
lmat), orbs%wfna(mbas, mo, 0:
lmat), orbs%wfnb(mbas, mo, 0:
lmat))
1091 ALLOCATE (orbs%pmat(mbas, mbas, 0:
lmat), orbs%pmata(mbas, mbas, 0:
lmat), orbs%pmatb(mbas, mbas, 0:
lmat))
1096 ALLOCATE (orbs%ener(mo, 0:
lmat), orbs%enera(mo, 0:
lmat), orbs%enerb(mo, 0:
lmat))
1101 ALLOCATE (orbs%refene(mo, 0:
lmat, 2), orbs%refchg(mo, 0:
lmat, 2), orbs%refnod(mo, 0:
lmat, 2))
1105 ALLOCATE (orbs%wrefene(mo, 0:
lmat, 2), orbs%wrefchg(mo, 0:
lmat, 2), orbs%wrefnod(mo, 0:
lmat, 2))
1106 orbs%wrefene = 0._dp
1107 orbs%wrefchg = 0._dp
1108 orbs%wrefnod = 0._dp
1109 ALLOCATE (orbs%crefene(mo, 0:
lmat, 2), orbs%crefchg(mo, 0:
lmat, 2), orbs%crefnod(mo, 0:
lmat, 2))
1110 orbs%crefene = 0._dp
1111 orbs%crefchg = 0._dp
1112 orbs%crefnod = 0._dp
1113 ALLOCATE (orbs%rcmax(mo, 0:
lmat, 2))
1115 ALLOCATE (orbs%wpsir0(mo, 2), orbs%tpsir0(mo, 2))
1118 ALLOCATE (orbs%reftype(mo, 0:
lmat, 2))
1130 cpassert(
ASSOCIATED(orbs))
1132 IF (
ASSOCIATED(orbs%wfn))
THEN
1133 DEALLOCATE (orbs%wfn, orbs%wfna, orbs%wfnb)
1135 IF (
ASSOCIATED(orbs%pmat))
THEN
1136 DEALLOCATE (orbs%pmat, orbs%pmata, orbs%pmatb)
1138 IF (
ASSOCIATED(orbs%ener))
THEN
1139 DEALLOCATE (orbs%ener, orbs%enera, orbs%enerb)
1141 IF (
ASSOCIATED(orbs%refene))
THEN
1142 DEALLOCATE (orbs%refene)
1144 IF (
ASSOCIATED(orbs%refchg))
THEN
1145 DEALLOCATE (orbs%refchg)
1147 IF (
ASSOCIATED(orbs%refnod))
THEN
1148 DEALLOCATE (orbs%refnod)
1150 IF (
ASSOCIATED(orbs%wrefene))
THEN
1151 DEALLOCATE (orbs%wrefene)
1153 IF (
ASSOCIATED(orbs%wrefchg))
THEN
1154 DEALLOCATE (orbs%wrefchg)
1156 IF (
ASSOCIATED(orbs%wrefnod))
THEN
1157 DEALLOCATE (orbs%wrefnod)
1159 IF (
ASSOCIATED(orbs%crefene))
THEN
1160 DEALLOCATE (orbs%crefene)
1162 IF (
ASSOCIATED(orbs%crefchg))
THEN
1163 DEALLOCATE (orbs%crefchg)
1165 IF (
ASSOCIATED(orbs%crefnod))
THEN
1166 DEALLOCATE (orbs%crefnod)
1168 IF (
ASSOCIATED(orbs%rcmax))
THEN
1169 DEALLOCATE (orbs%rcmax)
1171 IF (
ASSOCIATED(orbs%wpsir0))
THEN
1172 DEALLOCATE (orbs%wpsir0)
1174 IF (
ASSOCIATED(orbs%tpsir0))
THEN
1175 DEALLOCATE (orbs%tpsir0)
1177 IF (
ASSOCIATED(orbs%reftype))
THEN
1178 DEALLOCATE (orbs%reftype)
1194 REAL(kind=
dp),
INTENT(OUT) :: hf_frac
1195 LOGICAL,
INTENT(OUT) :: do_hfx
1198 INTEGER,
INTENT(IN) :: extype
1200 INTEGER :: i, j, nr, nu, pot_type
1201 REAL(kind=
dp) :: scale_coulomb, scale_longrange
1202 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: abscissa, weights
1206 IF (
ASSOCIATED(
atom%xc_section))
THEN
1207 xc_section =>
atom%xc_section
1212 atom%hfx_pot%scale_longrange = 0.0_dp
1213 atom%hfx_pot%scale_coulomb = 1.0_dp
1226 SELECT CASE (pot_type)
1228 cpwarn(
"Potential not implemented, use Coulomb instead!")
1230 atom%hfx_pot%scale_longrange = 0.0_dp
1231 atom%hfx_pot%scale_coulomb = scale_coulomb
1233 atom%hfx_pot%scale_coulomb = 0.0_dp
1234 atom%hfx_pot%scale_longrange = scale_longrange
1236 atom%hfx_pot%scale_coulomb = 1.0_dp
1237 atom%hfx_pot%scale_longrange = -1.0_dp
1239 atom%hfx_pot%scale_coulomb = scale_coulomb
1240 atom%hfx_pot%scale_longrange = scale_longrange
1246 cpabort(
"Only numerical and semi-analytic lrHF exchange available!")
1249 IF (
atom%hfx_pot%scale_longrange /= 0.0_dp .AND. extype ==
do_numeric .AND. .NOT.
ALLOCATED(
atom%hfx_pot%kernel))
THEN
1252 IF (
atom%hfx_pot%do_gh)
THEN
1256 ALLOCATE (weights(
atom%hfx_pot%nr_gh), abscissa(
atom%hfx_pot%nr_gh))
1257 CALL get_gauss_hermite_weights(abscissa, weights,
atom%hfx_pot%nr_gh)
1259 nr =
atom%basis%grid%nr
1260 ALLOCATE (
atom%hfx_pot%kernel(nr,
atom%hfx_pot%nr_gh, 0:
atom%state%maxl_calc +
atom%state%maxl_occ))
1261 atom%hfx_pot%kernel = 0.0_dp
1262 DO nu = 0,
atom%state%maxl_calc +
atom%state%maxl_occ
1263 DO i = 1,
atom%hfx_pot%nr_gh
1266 *abscissa(i)*
atom%basis%grid%rad(j), nu)*sqrt(weights(i))
1274 nr =
atom%basis%grid%nr
1275 ALLOCATE (
atom%hfx_pot%kernel(nr, nr, 0:
atom%state%maxl_calc +
atom%state%maxl_occ))
1276 atom%hfx_pot%kernel = 0.0_dp
1277 DO nu = 0,
atom%state%maxl_calc +
atom%state%maxl_occ
1281 *
atom%basis%grid%rad(i)*
atom%basis%grid%rad(j), nu)
1288 NULLIFY (xc_section)
1300 SUBROUTINE get_gauss_hermite_weights(abscissa, weights, nn)
1301 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: abscissa, weights
1302 INTEGER,
INTENT(IN) :: nn
1304 INTEGER :: counter, ii, info, liwork, lwork
1305 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: iwork
1306 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: diag, subdiag, work
1307 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: eigenvec
1311 ALLOCATE (work(1), iwork(1), diag(2*nn), subdiag(2*nn - 1), eigenvec(2*nn, 2*nn))
1316 subdiag(ii) = sqrt(real(ii, kind=
dp)/2.0_dp)
1320 CALL dstevd(
'V', 2*nn, diag, subdiag, eigenvec, 2*nn, work, lwork, iwork, liwork, info)
1323 cpabort(
'Finding size of working matrices failed!')
1327 lwork = int(work(1))
1329 DEALLOCATE (work, iwork)
1330 ALLOCATE (work(lwork), iwork(liwork))
1333 CALL dstevd(
'V', 2*nn, diag, subdiag, eigenvec, 2*nn, work, lwork, iwork, liwork, info)
1336 cpabort(
'Eigenvalue decomposition failed!')
1339 DEALLOCATE (work, iwork, subdiag)
1345 IF (diag(ii) > 0.0_dp)
THEN
1346 counter = counter + 1
1347 abscissa(counter) = diag(ii)
1348 weights(counter) =
rootpi*eigenvec(1, ii)**2
1351 IF (counter /= nn)
THEN
1352 cpabort(
'Have not found enough or too many zeros!')
1355 END SUBROUTINE get_gauss_hermite_weights
1365 INTEGER,
DIMENSION(0:lmat),
INTENT(IN) :: n
1366 INTEGER,
INTENT(IN),
OPTIONAL :: lmax
1371 IF (
PRESENT(lmax))
THEN
1377 cpassert(.NOT.
ASSOCIATED(opmat))
1382 ALLOCATE (opmat%op(m, m, 0:lm))
1394 cpassert(
ASSOCIATED(opmat))
1397 DEALLOCATE (opmat%op)
1414 cpassert(.NOT.
ASSOCIATED(opgrid))
1422 ALLOCATE (opgrid%op(nr))
1434 cpassert(
ASSOCIATED(opgrid))
1436 NULLIFY (opgrid%grid)
1437 DEALLOCATE (opgrid%op)
1453 INTEGER,
INTENT(IN) :: zval
1454 REAL(
dp),
INTENT(OUT) :: cval, aval
1455 INTEGER,
DIMENSION(0:lmat),
INTENT(OUT) :: ngto, ival
1468 cval = 2.14774520_dp
1469 aval = 0.04850670_dp
1472 cval = 2.08932430_dp
1473 aval = 0.02031060_dp
1476 cval = 2.09753060_dp
1477 aval = 0.03207070_dp
1480 cval = 2.10343410_dp
1481 aval = 0.03591970_dp
1485 cval = 2.10662820_dp
1486 aval = 0.05292410_dp
1490 cval = 2.13743840_dp
1491 aval = 0.06291970_dp
1495 cval = 2.08687310_dp
1496 aval = 0.08350860_dp
1500 cval = 2.12318180_dp
1501 aval = 0.09899170_dp
1505 cval = 2.13164810_dp
1506 aval = 0.11485350_dp
1510 cval = 2.11413310_dp
1511 aval = 0.00922630_dp
1516 cval = 2.12183620_dp
1517 aval = 0.01215850_dp
1522 cval = 2.06073230_dp
1523 aval = 0.01449350_dp
1528 cval = 2.08563660_dp
1529 aval = 0.01861460_dp
1534 cval = 2.04879270_dp
1535 aval = 0.02147790_dp
1540 cval = 2.06216660_dp
1541 aval = 0.01978920_dp
1546 cval = 2.04628670_dp
1547 aval = 0.02451470_dp
1552 cval = 2.08675200_dp
1553 aval = 0.02635040_dp
1558 cval = 2.02715220_dp
1559 aval = 0.01822040_dp
1564 cval = 2.01465650_dp
1565 aval = 0.01646570_dp
1570 cval = 2.01605240_dp
1571 aval = 0.01254190_dp
1577 cval = 2.01800000_dp
1578 aval = 0.01195490_dp
1585 cval = 1.98803560_dp
1586 aval = 0.02492140_dp
1593 cval = 1.98984000_dp
1594 aval = 0.02568400_dp
1601 cval = 2.01694380_dp
1602 aval = 0.02664480_dp
1609 cval = 2.01824090_dp
1610 aval = 0.01355000_dp
1617 cval = 1.98359400_dp
1618 aval = 0.01702210_dp
1625 cval = 1.96797340_dp
1626 aval = 0.02163180_dp
1633 cval = 1.98955180_dp
1634 aval = 0.02304480_dp
1641 cval = 1.98074320_dp
1642 aval = 0.02754320_dp
1649 cval = 2.00551070_dp
1650 aval = 0.02005530_dp
1657 cval = 2.00000030_dp
1658 aval = 0.02003000_dp
1665 cval = 2.00609100_dp
1666 aval = 0.02055620_dp
1673 cval = 2.00701000_dp
1674 aval = 0.02230400_dp
1681 cval = 2.01508710_dp
1682 aval = 0.02685790_dp
1689 cval = 2.01960430_dp
1690 aval = 0.02960430_dp
1697 cval = 2.00031000_dp
1698 aval = 0.00768400_dp
1706 cval = 1.99563960_dp
1707 aval = 0.01401940_dp
1714 cval = 1.98971210_dp
1715 aval = 0.01558470_dp
1721 cval = 1.97976190_dp
1722 aval = 0.01705520_dp
1728 cval = 1.97989290_dp
1729 aval = 0.01527040_dp
1735 cval = 1.97909240_dp
1736 aval = 0.01879720_dp
1742 cval = 1.98508430_dp
1743 aval = 0.01497550_dp
1750 cval = 1.98515010_dp
1751 aval = 0.01856670_dp
1758 cval = 1.98502970_dp
1759 aval = 0.01487000_dp
1766 cval = 1.97672850_dp
1767 aval = 0.01762500_dp
1775 cval = 1.97862730_dp
1776 aval = 0.01863310_dp
1783 cval = 1.97990020_dp
1784 aval = 0.01347150_dp
1791 cval = 1.97979410_dp
1792 aval = 0.00890265_dp
1799 cval = 1.98001000_dp
1800 aval = 0.00895215_dp
1807 cval = 1.97979980_dp
1808 aval = 0.01490290_dp
1815 cval = 1.98009310_dp
1816 aval = 0.01490390_dp
1823 cval = 1.97794750_dp
1824 aval = 0.01425880_dp
1832 cval = 1.97784450_dp
1833 aval = 0.01430130_dp
1841 cval = 1.97784450_dp
1842 aval = 0.00499318_dp
1850 cval = 1.97764820_dp
1851 aval = 0.00500392_dp
1859 cval = 1.97765150_dp
1860 aval = 0.00557083_dp
1868 cval = 1.97768750_dp
1869 aval = 0.00547531_dp
1879 cval = 1.96986600_dp
1880 aval = 0.00813143_dp
1890 cval = 1.97765720_dp
1891 aval = 0.00489201_dp
1901 cval = 1.97768120_dp
1902 aval = 0.00499000_dp
1912 cval = 1.97745700_dp
1913 aval = 0.00615587_dp
1923 cval = 1.97570240_dp
1924 aval = 0.00769959_dp
1934 cval = 1.97629350_dp
1935 aval = 0.00706610_dp
1945 cval = 1.96900000_dp
1946 aval = 0.01019150_dp
1956 cval = 1.97350000_dp
1957 aval = 0.01334320_dp
1967 cval = 1.97493000_dp
1968 aval = 0.01331360_dp
1977 cval = 1.97597670_dp
1978 aval = 0.01434040_dp
1988 cval = 1.97809240_dp
1989 aval = 0.01529430_dp
1999 cval = 1.97644360_dp
2000 aval = 0.01312770_dp
2010 cval = 1.96998000_dp
2011 aval = 0.01745150_dp
2021 cval = 1.97223830_dp
2022 aval = 0.01639750_dp
2032 cval = 1.97462110_dp
2033 aval = 0.01603680_dp
2043 cval = 1.97756000_dp
2044 aval = 0.02030570_dp
2054 cval = 1.97645760_dp
2055 aval = 0.02057180_dp
2065 cval = 1.97725820_dp
2066 aval = 0.02058210_dp
2076 cval = 1.97749380_dp
2077 aval = 0.02219380_dp
2087 cval = 1.97946280_dp
2088 aval = 0.02216280_dp
2098 cval = 1.97852130_dp
2099 aval = 0.02168500_dp
2109 cval = 1.98045190_dp
2110 aval = 0.02177860_dp
2120 cval = 1.97000000_dp
2121 aval = 0.02275000_dp
2131 cval = 1.97713580_dp
2132 aval = 0.02317030_dp
2142 cval = 1.97537880_dp
2143 aval = 0.02672860_dp
2153 cval = 1.97545360_dp
2154 aval = 0.02745360_dp
2164 cval = 1.97338370_dp
2165 aval = 0.02616310_dp
2175 cval = 1.97294240_dp
2176 aval = 0.02429220_dp
2186 cval = 1.98000000_dp
2187 aval = 0.01400000_dp
2197 cpabort(
"No geometrical basis set data are available for the selected atom number.")
2210 SUBROUTINE read_basis_set(element_symbol, basis, basis_set_name, basis_set_file, &
2213 CHARACTER(LEN=*),
INTENT(IN) :: element_symbol
2215 CHARACTER(LEN=*),
INTENT(IN) :: basis_set_name, basis_set_file
2219 INTEGER :: lmin, lmax
2220 INTEGER,
ALLOCATABLE :: l(:)
2221 REAL(
dp),
ALLOCATABLE :: zet(:), gcc(:, :)
2222 END TYPE basis_set_type
2224 CHARACTER(len=20*default_string_length) :: line_att
2225 CHARACTER(LEN=240) :: line
2226 CHARACTER(LEN=242) :: line2
2227 CHARACTER(LEN=LEN(basis_set_name)) :: bsname
2228 CHARACTER(LEN=LEN(basis_set_name) + 2) :: bsname2
2229 CHARACTER(LEN=LEN(element_symbol)) :: symbol
2230 CHARACTER(LEN=LEN(element_symbol) + 2) :: symbol2
2231 INTEGER :: i, ii, ipgf, iset, ishell, j, k, &
2232 lshell, n, nj, npgf, ns, nset, nshell, &
2234 INTEGER,
ALLOCATABLE :: shells(:)
2235 LOGICAL :: found, is_ok, match, read_from_input
2236 REAL(
dp) :: expzet, gcca, prefac, zeta
2237 TYPE(basis_set_type),
ALLOCATABLE :: sets(:)
2242 bsname = basis_set_name
2243 symbol = element_symbol
2245 read_from_input = .false.
2247 IF (read_from_input)
THEN
2254 CALL val_get(val, c_val=line_att)
2255 READ (line_att, *) nset
2268 line2 =
" "//line//
" "
2269 symbol2 =
" "//trim(symbol)//
" "
2270 bsname2 =
" "//trim(bsname)//
" "
2271 strlen1 = len_trim(symbol2) + 1
2272 strlen2 = len_trim(bsname2) + 1
2274 IF ((index(line2, symbol2(:strlen1)) > 0) .AND. &
2275 (index(line2, bsname2(:strlen2)) > 0)) match = .true.
2285 cpabort(
"End of file reached and the requested basis set was not found.")
2292 ALLOCATE (sets(nset))
2294 associate(set => sets(iset))
2295 IF (read_from_input)
THEN
2298 CALL val_get(val, c_val=line_att)
2299 READ (line_att, *) n, set%lmin, set%lmax, npgf
2309 cpassert(set%lmin >= 0 .AND. set%lmax >= set%lmin .AND. npgf > 0)
2310 ALLOCATE (shells(set%lmin:set%lmax))
2311 DO lshell = set%lmin, set%lmax
2312 IF (read_from_input)
THEN
2313 READ (line_att, *) shells(lshell)
2319 IF (read_from_input)
THEN
2320 cpassert(len_trim(line_att) == 0)
2322 cpassert(all(shells >= 0))
2323 nshell = sum(shells)
2324 ALLOCATE (set%l(nshell), set%zet(npgf), set%gcc(npgf, nshell))
2326 DO lshell = set%lmin, set%lmax
2327 set%l(ishell + 1:ishell + shells(lshell)) = lshell
2328 ishell = ishell + shells(lshell)
2332 IF (read_from_input)
THEN
2335 CALL val_get(val, c_val=line_att)
2336 READ (line_att, *) set%zet(ipgf), set%gcc(ipgf, :)
2339 DO ishell = 1, nshell
2352 DO j = sets(i)%lmin, min(sets(i)%lmax,
lmat)
2353 basis%nprim(j) = basis%nprim(j) +
SIZE(sets(i)%zet)
2355 DO j = 1,
SIZE(sets(i)%l)
2357 IF (k <=
lmat) basis%nbas(k) = basis%nbas(k) + 1
2361 nj = maxval(basis%nprim)
2362 ns = maxval(basis%nbas)
2363 ALLOCATE (basis%am(nj, 0:
lmat))
2365 ALLOCATE (basis%cm(nj, ns, 0:
lmat))
2372 IF (j >= sets(i)%lmin .AND. j <= sets(i)%lmax)
THEN
2373 DO ipgf = 1,
SIZE(sets(i)%zet)
2374 basis%am(nj + ipgf, j) = sets(i)%zet(ipgf)
2376 DO ii = 1,
SIZE(sets(i)%l)
2377 IF (sets(i)%l(ii) == j)
THEN
2379 DO ipgf = 1,
SIZE(sets(i)%zet)
2380 basis%cm(nj + ipgf, ns, j) = sets(i)%gcc(ipgf, ii)
2384 nj = nj +
SIZE(sets(i)%zet)
2391 expzet = 0.25_dp*real(2*j + 3,
dp)
2392 prefac = sqrt(
rootpi/2._dp**(j + 2)*
dfac(2*j + 1))
2393 DO ipgf = 1, basis%nprim(j)
2394 DO ii = 1, basis%nbas(j)
2395 gcca = basis%cm(ipgf, ii, j)
2396 zeta = 2._dp*basis%am(ipgf, j)
2397 basis%cm(ipgf, ii, j) = zeta**expzet*gcca/prefac
2402 END SUBROUTINE read_basis_set
2411 TYPE(section_vals_type),
POINTER :: opt_section
2413 INTEGER :: miter, ndiis
2414 REAL(kind=dp) :: damp, eps_diis, eps_scf
2416 CALL section_vals_val_get(opt_section,
"MAX_ITER", i_val=miter)
2417 CALL section_vals_val_get(opt_section,
"EPS_SCF", r_val=eps_scf)
2418 CALL section_vals_val_get(opt_section,
"N_DIIS", i_val=ndiis)
2419 CALL section_vals_val_get(opt_section,
"EPS_DIIS", r_val=eps_diis)
2420 CALL section_vals_val_get(opt_section,
"DAMPING", r_val=damp)
2422 optimization%max_iter = miter
2423 optimization%eps_scf = eps_scf
2424 optimization%n_diis = ndiis
2425 optimization%eps_diis = eps_diis
2426 optimization%damping = damp
2437 TYPE(section_vals_type),
POINTER :: potential_section
2438 INTEGER,
INTENT(IN) :: zval
2440 CHARACTER(LEN=default_path_length) :: pseudo_fn
2441 CHARACTER(LEN=default_string_length) :: pseudo_name
2443 REAL(dp),
DIMENSION(:),
POINTER :: convals
2444 TYPE(section_vals_type),
POINTER :: ecp_potential_section, &
2445 gth_potential_section
2448 CALL section_vals_val_get(potential_section,
"PSEUDO_TYPE", i_val=potential%ppot_type)
2450 SELECT CASE (potential%ppot_type)
2452 CALL section_vals_val_get(potential_section,
"POTENTIAL_FILE_NAME", c_val=pseudo_fn)
2453 CALL section_vals_val_get(potential_section,
"POTENTIAL_NAME", c_val=pseudo_name)
2454 gth_potential_section => section_vals_get_subs_vals(potential_section,
"GTH_POTENTIAL")
2455 CALL read_gth_potential(ptable(zval)%symbol, potential%gth_pot, &
2456 pseudo_name, pseudo_fn, gth_potential_section)
2458 CALL section_vals_val_get(potential_section,
"POTENTIAL_FILE_NAME", c_val=pseudo_fn)
2459 CALL section_vals_val_get(potential_section,
"POTENTIAL_NAME", c_val=pseudo_name)
2460 ecp_potential_section => section_vals_get_subs_vals(potential_section,
"ECP")
2462 pseudo_name, pseudo_fn, ecp_potential_section)
2464 CALL section_vals_val_get(potential_section,
"POTENTIAL_FILE_NAME", c_val=pseudo_fn)
2465 CALL section_vals_val_get(potential_section,
"POTENTIAL_NAME", c_val=pseudo_name)
2466 CALL atom_read_upf(potential%upf_pot, pseudo_fn)
2467 potential%upf_pot%pname = pseudo_name
2469 cpabort(
"Pseudopotential type SGP is not implemented.")
2473 cpabort(
"Invalid pseudopotential type selected. Check the code!")
2476 potential%ppot_type = no_pseudo
2481 CALL section_vals_val_get(potential_section,
"CONFINEMENT_TYPE", i_val=ic)
2482 potential%conf_type = ic
2483 IF (potential%conf_type == no_conf)
THEN
2484 potential%acon = 0.0_dp
2485 potential%rcon = 4.0_dp
2486 potential%scon = 2.0_dp
2487 potential%confinement = .false.
2488 ELSE IF (potential%conf_type == poly_conf)
THEN
2489 CALL section_vals_val_get(potential_section,
"CONFINEMENT", r_vals=convals)
2490 IF (
SIZE(convals) >= 1)
THEN
2491 IF (convals(1) > 0.0_dp)
THEN
2492 potential%confinement = .true.
2493 potential%acon = convals(1)
2494 IF (
SIZE(convals) >= 2)
THEN
2495 potential%rcon = convals(2)
2497 potential%rcon = 4.0_dp
2499 IF (
SIZE(convals) >= 3)
THEN
2500 potential%scon = convals(3)
2502 potential%scon = 2.0_dp
2505 potential%confinement = .false.
2508 potential%confinement = .false.
2510 ELSE IF (potential%conf_type == barrier_conf)
THEN
2511 potential%acon = 200.0_dp
2512 potential%rcon = 4.0_dp
2513 potential%scon = 12.0_dp
2514 potential%confinement = .true.
2515 CALL section_vals_val_get(potential_section,
"CONFINEMENT", r_vals=convals)
2516 IF (
SIZE(convals) >= 1)
THEN
2517 IF (convals(1) > 0.0_dp)
THEN
2518 potential%acon = convals(1)
2519 IF (
SIZE(convals) >= 2)
THEN
2520 potential%rcon = convals(2)
2522 IF (
SIZE(convals) >= 3)
THEN
2523 potential%scon = convals(3)
2526 potential%confinement = .false.
2539 potential%confinement = .false.
2541 CALL atom_release_upf(potential%upf_pot)
2552 SUBROUTINE read_gth_potential(element_symbol, potential, pseudo_name, pseudo_file, &
2555 CHARACTER(LEN=*),
INTENT(IN) :: element_symbol
2557 CHARACTER(LEN=*),
INTENT(IN) :: pseudo_name, pseudo_file
2558 TYPE(section_vals_type),
POINTER :: potential_section
2560 CHARACTER(LEN=240) :: line
2561 CHARACTER(LEN=242) :: line2
2562 CHARACTER(len=5*default_string_length) :: line_att
2563 CHARACTER(LEN=LEN(element_symbol)) :: symbol
2564 CHARACTER(LEN=LEN(element_symbol)+2) :: symbol2
2565 CHARACTER(LEN=LEN(pseudo_name)) :: apname
2566 CHARACTER(LEN=LEN(pseudo_name)+2) :: apname2
2567 INTEGER :: i, ic, ipot, j, l, nlmax, strlen1, &
2569 INTEGER,
DIMENSION(0:lmat) :: elec_conf
2570 LOGICAL :: found, is_ok, match, read_from_input
2571 TYPE(cp_sll_val_type),
POINTER ::
list
2572 TYPE(val_type),
POINTER :: val
2576 apname = pseudo_name
2577 symbol = element_symbol
2579 potential%symbol = symbol
2580 potential%pname = apname
2582 potential%rc = 0._dp
2584 potential%cl = 0._dp
2586 potential%rcnl = 0._dp
2587 potential%hnl = 0._dp
2588 potential%soc = .false.
2589 potential%knl = 0._dp
2591 potential%lpotextended = .false.
2592 potential%lsdpot = .false.
2593 potential%nlcc = .false.
2594 potential%nexp_lpot = 0
2595 potential%nexp_lsd = 0
2596 potential%nexp_nlcc = 0
2598 read_from_input = .false.
2599 CALL section_vals_get(potential_section, explicit=read_from_input)
2600 IF (read_from_input)
THEN
2601 CALL section_vals_list_get(potential_section,
"_DEFAULT_KEYWORD_",
list=
list)
2602 CALL uppercase(symbol)
2603 CALL uppercase(apname)
2606 is_ok = cp_sll_val_next(
list, val)
2608 CALL val_get(val, c_val=line_att)
2609 READ (line_att, *) elec_conf(l)
2610 CALL remove_word(line_att)
2611 DO WHILE (len_trim(line_att) /= 0)
2613 READ (line_att, *) elec_conf(l)
2614 CALL remove_word(line_att)
2616 potential%econf(0:
lmat) = elec_conf(0:
lmat)
2617 potential%zion = real(sum(elec_conf), dp)
2619 is_ok = cp_sll_val_next(
list, val)
2621 CALL val_get(val, c_val=line_att)
2622 READ (line_att, *) potential%rc
2623 CALL remove_word(line_att)
2625 READ (line_att, *) potential%ncl
2626 CALL remove_word(line_att)
2627 DO i = 1, potential%ncl
2628 READ (line_att, *) potential%cl(i)
2629 CALL remove_word(line_att)
2633 is_ok = cp_sll_val_next(
list, val)
2635 CALL val_get(val, c_val=line_att)
2636 IF (index(line_att,
"LPOT") /= 0)
THEN
2637 potential%lpotextended = .true.
2638 CALL remove_word(line_att)
2639 READ (line_att, *) potential%nexp_lpot
2640 DO ipot = 1, potential%nexp_lpot
2641 is_ok = cp_sll_val_next(
list, val)
2643 CALL val_get(val, c_val=line_att)
2644 READ (line_att, *) potential%alpha_lpot(ipot)
2645 CALL remove_word(line_att)
2646 READ (line_att, *) potential%nct_lpot(ipot)
2647 CALL remove_word(line_att)
2648 DO ic = 1, potential%nct_lpot(ipot)
2649 READ (line_att, *) potential%cval_lpot(ic, ipot)
2650 CALL remove_word(line_att)
2653 ELSE IF (index(line_att,
"NLCC") /= 0)
THEN
2654 potential%nlcc = .true.
2655 CALL remove_word(line_att)
2656 READ (line_att, *) potential%nexp_nlcc
2657 DO ipot = 1, potential%nexp_nlcc
2658 is_ok = cp_sll_val_next(
list, val)
2660 CALL val_get(val, c_val=line_att)
2661 READ (line_att, *) potential%alpha_nlcc(ipot)
2662 CALL remove_word(line_att)
2663 READ (line_att, *) potential%nct_nlcc(ipot)
2664 CALL remove_word(line_att)
2665 DO ic = 1, potential%nct_nlcc(ipot)
2666 READ (line_att, *) potential%cval_nlcc(ic, ipot)
2668 potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2669 CALL remove_word(line_att)
2672 ELSE IF (index(line_att,
"LSD") /= 0)
THEN
2673 potential%lsdpot = .true.
2674 CALL remove_word(line_att)
2675 READ (line_att, *) potential%nexp_lsd
2676 DO ipot = 1, potential%nexp_lsd
2677 is_ok = cp_sll_val_next(
list, val)
2679 CALL val_get(val, c_val=line_att)
2680 READ (line_att, *) potential%alpha_lsd(ipot)
2681 CALL remove_word(line_att)
2682 READ (line_att, *) potential%nct_lsd(ipot)
2683 CALL remove_word(line_att)
2684 DO ic = 1, potential%nct_lsd(ipot)
2685 READ (line_att, *) potential%cval_lsd(ic, ipot)
2686 CALL remove_word(line_att)
2694 READ (line_att, *) nlmax
2695 CALL remove_word(line_att)
2696 IF (index(line_att,
"SOC") /= 0) potential%soc = .true.
2700 is_ok = cp_sll_val_next(
list, val)
2702 CALL val_get(val, c_val=line_att)
2703 READ (line_att, *) potential%rcnl(l)
2704 CALL remove_word(line_att)
2705 READ (line_att, *) potential%nl(l)
2706 CALL remove_word(line_att)
2707 DO i = 1, potential%nl(l)
2709 READ (line_att, *) potential%hnl(1, 1, l)
2710 CALL remove_word(line_att)
2712 cpassert(len_trim(line_att) == 0)
2713 is_ok = cp_sll_val_next(
list, val)
2715 CALL val_get(val, c_val=line_att)
2716 READ (line_att, *) potential%hnl(i, i, l)
2717 CALL remove_word(line_att)
2719 DO j = i + 1, potential%nl(l)
2720 READ (line_att, *) potential%hnl(i, j, l)
2721 potential%hnl(j, i, l) = potential%hnl(i, j, l)
2722 CALL remove_word(line_att)
2725 IF (potential%soc .AND. l /= 0)
THEN
2726 is_ok = cp_sll_val_next(
list, val)
2728 CALL val_get(val, c_val=line_att)
2729 DO i = 1, potential%nl(l)
2731 READ (line_att, *) potential%knl(1, 1, l)
2732 CALL remove_word(line_att)
2734 cpassert(len_trim(line_att) == 0)
2735 is_ok = cp_sll_val_next(
list, val)
2737 CALL val_get(val, c_val=line_att)
2738 READ (line_att, *) potential%knl(i, i, l)
2739 CALL remove_word(line_att)
2741 DO j = i + 1, potential%nl(l)
2742 READ (line_att, *) potential%knl(i, j, l)
2743 potential%knl(j, i, l) = potential%knl(i, j, l)
2744 CALL remove_word(line_att)
2748 cpassert(len_trim(line_att) == 0)
2753 TYPE(cp_parser_type) :: parser
2754 CALL parser_create(parser, pseudo_file)
2757 CALL parser_search_string(parser, trim(apname), .true., found, line)
2759 CALL uppercase(symbol)
2760 CALL uppercase(apname)
2763 CALL uppercase(line)
2764 line2 =
" "//line//
" "
2765 symbol2 =
" "//trim(symbol)//
" "
2766 apname2 =
" "//trim(apname)//
" "
2767 strlen1 = len_trim(symbol2) + 1
2768 strlen2 = len_trim(apname2) + 1
2770 IF ((index(line2, symbol2(:strlen1)) > 0) .AND. &
2771 (index(line2, apname2(:strlen2)) > 0)) match = .true.
2776 CALL parser_get_object(parser, elec_conf(l), newline=.true.)
2777 DO WHILE (parser_test_next_token(parser) ==
"INT")
2779 CALL parser_get_object(parser, elec_conf(l))
2781 potential%econf(0:
lmat) = elec_conf(0:
lmat)
2782 potential%zion = real(sum(elec_conf), dp)
2784 CALL parser_get_object(parser, potential%rc, newline=.true.)
2786 CALL parser_get_object(parser, potential%ncl)
2787 DO i = 1, potential%ncl
2788 CALL parser_get_object(parser, potential%cl(i))
2792 CALL parser_get_next_line(parser, 1)
2793 IF (parser_test_next_token(parser) ==
"INT")
THEN
2795 ELSE IF (parser_test_next_token(parser) ==
"STR")
THEN
2796 CALL parser_get_object(parser, line)
2797 IF (index(line,
"LPOT") /= 0)
THEN
2799 potential%lpotextended = .true.
2800 CALL parser_get_object(parser, potential%nexp_lpot)
2801 DO ipot = 1, potential%nexp_lpot
2802 CALL parser_get_object(parser, potential%alpha_lpot(ipot), newline=.true.)
2803 CALL parser_get_object(parser, potential%nct_lpot(ipot))
2804 DO ic = 1, potential%nct_lpot(ipot)
2805 CALL parser_get_object(parser, potential%cval_lpot(ic, ipot))
2808 ELSE IF (index(line,
"NLCC") /= 0)
THEN
2810 potential%nlcc = .true.
2811 CALL parser_get_object(parser, potential%nexp_nlcc)
2812 DO ipot = 1, potential%nexp_nlcc
2813 CALL parser_get_object(parser, potential%alpha_nlcc(ipot), newline=.true.)
2814 CALL parser_get_object(parser, potential%nct_nlcc(ipot))
2815 DO ic = 1, potential%nct_nlcc(ipot)
2816 CALL parser_get_object(parser, potential%cval_nlcc(ic, ipot))
2818 potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2821 ELSE IF (index(line,
"LSD") /= 0)
THEN
2823 potential%lsdpot = .true.
2824 CALL parser_get_object(parser, potential%nexp_lsd)
2825 DO ipot = 1, potential%nexp_lsd
2826 CALL parser_get_object(parser, potential%alpha_lsd(ipot), newline=.true.)
2827 CALL parser_get_object(parser, potential%nct_lsd(ipot))
2828 DO ic = 1, potential%nct_lsd(ipot)
2829 CALL parser_get_object(parser, potential%cval_lsd(ic, ipot))
2833 cpabort(
"Parsing of extended potential type failed.")
2836 cpabort(
"Invalid input token found.")
2840 CALL parser_get_object(parser, nlmax)
2842 IF (parser_test_next_token(parser) ==
"STR")
THEN
2843 CALL parser_get_object(parser, line)
2844 IF (index(line,
"SOC") /= 0) potential%soc = .true.
2848 CALL parser_get_object(parser, potential%rcnl(l), newline=.true.)
2849 CALL parser_get_object(parser, potential%nl(l))
2850 DO i = 1, potential%nl(l)
2852 CALL parser_get_object(parser, potential%hnl(i, i, l))
2854 CALL parser_get_object(parser, potential%hnl(i, i, l), newline=.true.)
2856 DO j = i + 1, potential%nl(l)
2857 CALL parser_get_object(parser, potential%hnl(i, j, l))
2858 potential%hnl(j, i, l) = potential%hnl(i, j, l)
2861 IF (potential%soc .AND. l /= 0)
THEN
2862 DO i = 1, potential%nl(l)
2863 CALL parser_get_object(parser, potential%knl(i, i, l), newline=.true.)
2864 DO j = i + 1, potential%nl(l)
2865 CALL parser_get_object(parser, potential%knl(i, j, l))
2866 potential%knl(j, i, l) = potential%knl(i, j, l)
2876 cpabort(
"End of file reached unexpectedly")
2881 CALL parser_release(parser)
2885 END SUBROUTINE read_gth_potential
2895 SUBROUTINE read_ecp_potential_file(element_symbol, potential, pseudo_name, pseudo_file, &
2896 potential_section, potential_found)
2898 CHARACTER(LEN=*),
INTENT(IN) :: element_symbol
2900 CHARACTER(LEN=*),
INTENT(IN) :: pseudo_name, pseudo_file
2901 TYPE(section_vals_type),
POINTER :: potential_section
2902 LOGICAL,
INTENT(OUT),
OPTIONAL :: potential_found
2904 CHARACTER(LEN=240) :: line
2905 CHARACTER(len=5*default_string_length) :: line_att
2906 CHARACTER(LEN=LEN(element_symbol)+1) :: symbol
2907 CHARACTER(LEN=LEN(pseudo_name)) :: apname
2908 INTEGER :: i, ic, l, ncore, nel
2909 LOGICAL :: found, is_ok, read_from_input
2910 TYPE(cp_sll_val_type),
POINTER :: list
2911 TYPE(val_type),
POINTER :: val
2913 apname = pseudo_name
2914 symbol = element_symbol
2915 IF (
PRESENT(potential_found)) potential_found = .false.
2916 CALL get_ptable_info(symbol, number=ncore)
2918 potential%symbol = symbol
2919 potential%pname = apname
2925 potential%aloc = 0.0_dp
2926 potential%bloc = 0.0_dp
2929 potential%apot = 0.0_dp
2930 potential%bpot = 0.0_dp
2932 read_from_input = .false.
2933 CALL section_vals_get(potential_section, explicit=read_from_input)
2934 IF (read_from_input)
THEN
2935 CALL section_vals_list_get(potential_section,
"_DEFAULT_KEYWORD_",
list=
list)
2937 is_ok = cp_sll_val_next(
list, val)
2939 CALL val_get(val, c_val=line_att)
2940 CALL remove_word(line_att)
2941 CALL remove_word(line_att)
2943 READ (line_att, *) nel
2944 potential%zion = real(ncore - nel, kind=dp)
2946 is_ok = cp_sll_val_next(
list, val)
2948 CALL val_get(val, c_val=line_att)
2950 IF (.NOT. cp_sll_val_next(
list, val))
EXIT
2951 CALL val_get(val, c_val=line_att)
2952 IF (index(line_att, element_symbol) == 0)
THEN
2953 potential%nloc = potential%nloc + 1
2955 READ (line_att, *) potential%nrloc(ic), potential%bloc(ic), potential%aloc(ic)
2962 CALL val_get(val, c_val=line_att)
2963 IF (index(line_att, element_symbol) == 0)
THEN
2964 potential%npot(l) = potential%npot(l) + 1
2965 ic = potential%npot(l)
2966 READ (line_att, *) potential%nrpot(ic, l), potential%bpot(ic, l), potential%apot(ic, l)
2968 potential%lmax = potential%lmax + 1
2971 IF (.NOT. cp_sll_val_next(
list, val))
EXIT
2976 TYPE(cp_parser_type) :: parser
2977 CALL parser_create(parser, pseudo_file)
2980 CALL parser_search_string(parser, trim(apname), .true., found, line)
2983 CALL parser_get_object(parser, line, newline=.true.)
2984 IF (trim(line) == element_symbol)
THEN
2985 IF (
PRESENT(potential_found)) potential_found = .true.
2986 CALL parser_get_object(parser, line, lower_to_upper=.true.)
2987 cpassert(trim(line) ==
"NELEC")
2989 CALL parser_get_object(parser, nel)
2990 potential%zion = real(ncore - nel, kind=dp)
2992 CALL parser_get_object(parser, line, newline=.true.)
2995 CALL parser_read_line(parser, 1)
2996 IF (parser_test_next_token(parser) ==
"STR")
EXIT
2997 potential%nloc = potential%nloc + 1
2999 CALL parser_get_object(parser, potential%nrloc(ic))
3000 CALL parser_get_object(parser, potential%bloc(ic))
3001 CALL parser_get_object(parser, potential%aloc(ic))
3005 CALL parser_get_object(parser, symbol)
3006 IF (symbol == element_symbol)
THEN
3008 potential%lmax = potential%lmax + 1
3010 CALL parser_read_line(parser, 1)
3011 IF (parser_test_next_token(parser) ==
"STR")
EXIT
3012 potential%npot(l) = potential%npot(l) + 1
3013 ic = potential%npot(l)
3014 CALL parser_get_object(parser, potential%nrpot(ic, l))
3015 CALL parser_get_object(parser, potential%bpot(ic, l))
3016 CALL parser_get_object(parser, potential%apot(ic, l))
3023 ELSE IF (line ==
"END")
THEN
3024 IF (
PRESENT(potential_found))
THEN
3025 CALL parser_release(parser)
3028 cpabort(
"Element not found in ECP library")
3032 IF (
PRESENT(potential_found))
THEN
3033 CALL parser_release(parser)
3036 cpabort(
"ECP type not found in library")
3041 CALL parser_release(parser)
3045 IF (read_from_input .AND.
PRESENT(potential_found)) potential_found = .true.
3048 potential%econf(0:3) = ptable(ncore)%e_conv(0:3)
3051 cpabort(
"Unknown Core State")
3054 potential%econf(0:3) = potential%econf(0:3) - ptable(2)%e_conv(0:3)
3056 potential%econf(0:3) = potential%econf(0:3) - ptable(10)%e_conv(0:3)
3058 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
3060 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
3061 potential%econf(2) = potential%econf(2) - 10
3063 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3065 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3066 potential%econf(2) = potential%econf(2) - 10
3068 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3070 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3071 potential%econf(2) = potential%econf(2) - 10
3072 potential%econf(3) = potential%econf(3) - 14
3074 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3075 potential%econf(3) = potential%econf(3) - 14
3077 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3078 potential%econf(2) = potential%econf(2) - 10
3079 potential%econf(3) = potential%econf(3) - 14
3082 cpassert(all(potential%econf >= 0))
3084 END SUBROUTINE read_ecp_potential_file
3094 SUBROUTINE read_ecp_potential_files(element_symbol, potential, pseudo_name, pseudo_files, &
3097 CHARACTER(LEN=*),
INTENT(IN) :: element_symbol
3099 CHARACTER(LEN=*),
INTENT(IN) :: pseudo_name
3100 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(IN) :: pseudo_files
3101 TYPE(section_vals_type),
POINTER :: potential_section
3103 CHARACTER(LEN=:),
ALLOCATABLE :: file_list
3105 LOGICAL :: potential_found
3107 DO i = 1,
SIZE(pseudo_files)
3108 CALL read_ecp_potential_file(element_symbol, potential, pseudo_name, pseudo_files(i), &
3109 potential_section, potential_found)
3110 IF (potential_found)
RETURN
3114 DO i = 1,
SIZE(pseudo_files)
3115 file_list = trim(file_list)//
"<"//trim(pseudo_files(i))//
"> "
3117 CALL cp_abort(__location__, &
3118 "The requested ECP potential <"//trim(pseudo_name)// &
3119 "> for element <"//trim(element_symbol)// &
3120 "> was not found in the potential files "//trim(file_list))
3122 END SUBROUTINE read_ecp_potential_files
3130 TYPE(grid_atom_type) :: grid1, grid2
3134 REAL(kind=dp) :: dr, dw
3137 IF (grid1%nr == grid2%nr)
THEN
3139 dr = abs(grid1%rad(i) - grid2%rad(i))
3140 dw = abs(grid1%wr(i) - grid2%wr(i))
3141 IF (dr + dw > 1.0e-12_dp)
THEN
Define the atom type and its sub types.
integer, parameter, public num_basis
subroutine, public read_atom_opt_section(optimization, opt_section)
...
subroutine, public create_atom_type(atom)
...
integer, parameter, public cgto_basis
integer, parameter, public gto_basis
integer, parameter, public sto_basis
subroutine, public release_atom_type(atom)
...
subroutine, public release_opmat(opmat)
...
subroutine, public release_atom_potential(potential)
...
subroutine, public init_atom_basis(basis, basis_section, zval, btyp)
Initialize the basis for the atomic code.
logical function, public atom_compare_grids(grid1, grid2)
...
subroutine, public release_opgrid(opgrid)
...
subroutine, public release_atom_orbs(orbs)
...
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 init_atom_potential(potential, potential_section, zval)
...
subroutine, public release_atom_basis(basis)
...
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
subroutine, public init_atom_basis_default_pp(basis)
...
subroutine, public create_opmat(opmat, n, lmax)
...
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
subroutine, public setup_hf_section(hf_frac, do_hfx, atom, xc_section, extype)
...
subroutine, public clementi_geobas(zval, cval, aval, ngto, ival)
...
subroutine, public create_opgrid(opgrid, grid)
...
Routines that process Quantum Espresso UPF files.
pure subroutine, public atom_release_upf(upfpot)
...
subroutine, public atom_read_upf(pot, upf_filename, read_header, read_augmentation)
...
Calculates Bessel functions.
elemental impure real(kind=dp) function, public bessel0(x, l)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public limpanuparb2011
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_read_line(parser, nline, at_end)
Read the next line from a logical unit "unit" (I/O node only). Skip (nline-1) lines and skip also all...
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
character(len=3) function, public parser_test_next_token(parser, string_length)
Test next input object.
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), parameter, public rootpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
subroutine, public deallocate_grid_atom(grid_atom)
Deallocate a Gaussian-type orbital (GTO) basis set data set.
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)
...
Utilities for string manipulations.
character(len=1), parameter, public newline
subroutine, public remove_word(string)
remove a word from a string (words are separated by white spaces)
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Provides all information about a basis set.
Provides all information about a pseudopotential.
Provides info about hartree-fock exchange (For now, we only support potentials that can be represente...
Information on optimization procedure.
Holds atomic orbitals and energies.
Provides all information on states and occupation.
Provides all information about an atomic kind.