(git:07c9450)
Loading...
Searching...
No Matches
atom_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Define the atom type and its sub types
10!> \author jgh
11!> \date 03.03.2008
12!> \version 1.0
13!>
14! **************************************************************************************************
16 USE atom_upf, ONLY: atom_read_upf,&
19 USE bessel_lib, ONLY: bessel0
20 USE bibliography, ONLY: limpanuparb2011,&
21 cite_reference
32 USE input_constants, ONLY: &
42 USE input_val_types, ONLY: val_get,&
44 USE kinds, ONLY: default_string_length,&
45 dp
46 USE mathconstants, ONLY: dfac,&
47 fac,&
48 pi,&
49 rootpi
51 ptable
56 USE string_utilities, ONLY: remove_word,&
58#include "./base/base_uses.f90"
59
60 IMPLICIT NONE
61
62 PRIVATE
63
64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_types'
65
66 ! maximum l-quantum number considered in atomic code/basis
67 INTEGER, PARAMETER :: lmat = 5
68
69 INTEGER, PARAMETER :: gto_basis = 100, &
70 cgto_basis = 101, &
71 sto_basis = 102, &
72 num_basis = 103
73
74 INTEGER, PARAMETER :: nmax = 25
75
76!> \brief Provides all information about a basis set
77! **************************************************************************************************
79 INTEGER :: basis_type = gto_basis
80 INTEGER, DIMENSION(0:lmat) :: nbas = 0
81 INTEGER, DIMENSION(0:lmat) :: nprim = 0
82 REAL(kind=dp), DIMENSION(:, :), POINTER :: am => null() !GTO exponents
83 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: cm => null() !Contraction coeffs
84 REAL(kind=dp), DIMENSION(:, :), POINTER :: as => null() !STO exponents
85 INTEGER, DIMENSION(:, :), POINTER :: ns => null() !STO n-quantum numbers
86 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: bf => null() !num. bsf
87 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dbf => null() !derivatives (num)
88 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ddbf => null() !2nd derivatives (num)
89 REAL(kind=dp) :: eps_eig = 0.0_dp
90 TYPE(grid_atom_type), POINTER :: grid => null()
91 LOGICAL :: geometrical = .false.
92 REAL(kind=dp) :: aval = 0.0_dp, cval = 0.0_dp
93 INTEGER, DIMENSION(0:lmat) :: start = 0
94 END TYPE atom_basis_type
95
96!> \brief Provides all information about a pseudopotential
97! **************************************************************************************************
99 CHARACTER(LEN=2) :: symbol = ""
100 CHARACTER(LEN=default_string_length) :: pname = ""
101 INTEGER, DIMENSION(0:lmat) :: econf = 0
102 REAL(dp) :: zion = 0.0_dp
103 REAL(dp) :: rc = 0.0_dp
104 INTEGER :: ncl = 0
105 REAL(dp), DIMENSION(5) :: cl = 0.0_dp
106 INTEGER, DIMENSION(0:lmat) :: nl = 0
107 REAL(dp), DIMENSION(0:lmat) :: rcnl = 0.0_dp
108 REAL(dp), DIMENSION(4, 4, 0:lmat) :: hnl = 0.0_dp
109 ! SOC
110 LOGICAL :: soc = .false.
111 REAL(dp), DIMENSION(4, 4, 0:lmat) :: knl = 0.0_dp
112 ! type extensions
113 ! NLCC
114 LOGICAL :: nlcc = .false.
115 INTEGER :: nexp_nlcc = 0
116 REAL(kind=dp), DIMENSION(10) :: alpha_nlcc = 0.0_dp
117 INTEGER, DIMENSION(10) :: nct_nlcc = 0
118 REAL(kind=dp), DIMENSION(4, 10) :: cval_nlcc = 0.0_dp
119 ! LSD potential
120 LOGICAL :: lsdpot = .false.
121 INTEGER :: nexp_lsd = 0
122 REAL(kind=dp), DIMENSION(10) :: alpha_lsd = 0.0_dp
123 INTEGER, DIMENSION(10) :: nct_lsd = 0
124 REAL(kind=dp), DIMENSION(4, 10) :: cval_lsd = 0.0_dp
125 ! extended local potential
126 LOGICAL :: lpotextended = .false.
127 INTEGER :: nexp_lpot = 0
128 REAL(kind=dp), DIMENSION(10) :: alpha_lpot = 0.0_dp
129 INTEGER, DIMENSION(10) :: nct_lpot = 0
130 REAL(kind=dp), DIMENSION(4, 10) :: cval_lpot = 0.0_dp
131 END TYPE atom_gthpot_type
132
134 CHARACTER(LEN=2) :: symbol = ""
135 CHARACTER(LEN=default_string_length) :: pname = ""
136 INTEGER, DIMENSION(0:lmat) :: econf = 0
137 REAL(dp) :: zion = 0.0_dp
138 INTEGER :: lmax = 0
139 INTEGER :: nloc = 0 ! # terms
140 INTEGER, DIMENSION(1:15) :: nrloc = 0 ! r**(n-2)
141 REAL(dp), DIMENSION(1:15) :: aloc = 0.0_dp ! coefficient
142 REAL(dp), DIMENSION(1:15) :: bloc = 0.0_dp ! exponent
143 INTEGER, DIMENSION(0:10) :: npot = 0 ! # terms
144 INTEGER, DIMENSION(1:15, 0:10) :: nrpot = 0 ! r**(n-2)
145 REAL(dp), DIMENSION(1:15, 0:10) :: apot = 0.0_dp ! coefficient
146 REAL(dp), DIMENSION(1:15, 0:10) :: bpot = 0.0_dp ! exponent
147 END TYPE atom_ecppot_type
148
150 CHARACTER(LEN=2) :: symbol = ""
151 CHARACTER(LEN=default_string_length) :: pname = ""
152 INTEGER, DIMENSION(0:lmat) :: econf = 0
153 REAL(dp) :: zion = 0.0_dp
154 INTEGER :: lmax = 0
155 LOGICAL :: has_nonlocal = .false.
156 INTEGER :: n_nonlocal = 0
157 LOGICAL, DIMENSION(0:5) :: is_nonlocal = .false.
158 REAL(kind=dp), DIMENSION(nmax) :: a_nonlocal = 0.0_dp
159 REAL(kind=dp), DIMENSION(nmax, 0:lmat) :: h_nonlocal = 0.0_dp
160 REAL(kind=dp), DIMENSION(nmax, nmax, 0:lmat) :: c_nonlocal = 0.0_dp
161 INTEGER :: n_local = 0
162 REAL(kind=dp) :: ac_local = 0.0_dp
163 REAL(kind=dp), DIMENSION(nmax) :: a_local = 0.0_dp
164 REAL(kind=dp), DIMENSION(nmax) :: c_local = 0.0_dp
165 LOGICAL :: has_nlcc = .false.
166 INTEGER :: n_nlcc = 0
167 REAL(kind=dp), DIMENSION(nmax) :: a_nlcc = 0.0_dp
168 REAL(kind=dp), DIMENSION(nmax) :: c_nlcc = 0.0_dp
169 END TYPE atom_sgppot_type
170
172 INTEGER :: ppot_type = 0
173 LOGICAL :: confinement = .false.
174 INTEGER :: conf_type = 0
175 REAL(dp) :: acon = 0.0_dp
176 REAL(dp) :: rcon = 0.0_dp
177 REAL(dp) :: scon = 0.0_dp
182 END TYPE atom_potential_type
183
184!> \brief Provides info about hartree-fock exchange (For now, we only support potentials that can be represented
185!> with Coulomb and longrange-coulomb potential)
186! **************************************************************************************************
188 REAL(kind=dp) :: scale_coulomb = 0.0_dp
189 REAL(kind=dp) :: scale_longrange = 0.0_dp
190 REAL(kind=dp) :: omega = 0.0_dp
191 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kernel
192 LOGICAL :: do_gh = .false.
193 INTEGER :: nr_gh = 0
194 END TYPE atom_hfx_type
195
196!> \brief Provides all information on states and occupation
197! **************************************************************************************************
199 REAL(kind=dp), DIMENSION(0:lmat, 10) :: occ = 0.0_dp
200 REAL(kind=dp), DIMENSION(0:lmat, 10) :: core = 0.0_dp
201 REAL(kind=dp), DIMENSION(0:lmat, 10) :: occupation = 0.0_dp
202 INTEGER :: maxl_occ = 0
203 INTEGER, DIMENSION(0:lmat) :: maxn_occ = 0
204 INTEGER :: maxl_calc = 0
205 INTEGER, DIMENSION(0:lmat) :: maxn_calc = 0
206 INTEGER :: multiplicity = 0
207 REAL(kind=dp), DIMENSION(0:lmat, 10) :: occa = 0.0_dp, occb = 0.0_dp
208 END TYPE atom_state
209
210!> \brief Holds atomic integrals
211! **************************************************************************************************
212 TYPE eri
213 REAL(kind=dp), DIMENSION(:, :), POINTER :: int => null()
214 END TYPE eri
215
217 INTEGER :: status = 0
218 INTEGER :: ppstat = 0
219 LOGICAL :: eri_coulomb = .false.
220 LOGICAL :: eri_exchange = .false.
221 LOGICAL :: all_nu = .false.
222 INTEGER, DIMENSION(0:lmat) :: n = 0, nne = 0
223 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ovlp => null(), kin => null(), core => null(), clsd => null()
224 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: utrans => null(), uptrans => null()
225 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: hnl => null()
226 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: conf => null()
227 TYPE(eri), DIMENSION(100) :: ceri = eri()
228 TYPE(eri), DIMENSION(100) :: eeri = eri()
229 INTEGER :: dkhstat = 0
230 INTEGER :: zorastat = 0
231 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: tzora => null()
232 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: hdkh => null()
233 END TYPE atom_integrals
234
235!> \brief Holds atomic orbitals and energies
236! **************************************************************************************************
238 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: wfn => null(), wfna => null(), wfnb => null()
239 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pmat => null(), pmata => null(), pmatb => null()
240 REAL(kind=dp), DIMENSION(:, :), POINTER :: ener => null(), enera => null(), enerb => null()
241 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: refene => null(), refchg => null(), refnod => null()
242 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: wrefene => null(), wrefchg => null(), wrefnod => null()
243 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: crefene => null(), crefchg => null(), crefnod => null()
244 REAL(kind=dp), DIMENSION(:, :), POINTER :: wpsir0 => null(), tpsir0 => null()
245 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rcmax => null()
246 CHARACTER(LEN=2), DIMENSION(:, :, :), POINTER :: reftype => null()
247 END TYPE atom_orbitals
248
249!> \brief Operator matrices
250! **************************************************************************************************
252 INTEGER, DIMENSION(0:lmat) :: n = 0
253 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: op => null()
254 END TYPE opmat_type
255
256!> \brief Operator grids
257! **************************************************************************************************
259 REAL(kind=dp), DIMENSION(:), POINTER :: op => null()
260 TYPE(grid_atom_type), POINTER :: grid => null()
261 END TYPE opgrid_type
262
263!> \brief All energies
264! **************************************************************************************************
265 TYPE atom_energy_type
266 REAL(kind=dp) :: etot = 0.0_dp
267 REAL(kind=dp) :: eband = 0.0_dp
268 REAL(kind=dp) :: ekin = 0.0_dp
269 REAL(kind=dp) :: epot = 0.0_dp
270 REAL(kind=dp) :: ecore = 0.0_dp
271 REAL(kind=dp) :: elsd = 0.0_dp
272 REAL(kind=dp) :: epseudo = 0.0_dp
273 REAL(kind=dp) :: eploc = 0.0_dp
274 REAL(kind=dp) :: epnl = 0.0_dp
275 REAL(kind=dp) :: exc = 0.0_dp
276 REAL(kind=dp) :: ecoulomb = 0.0_dp
277 REAL(kind=dp) :: eexchange = 0.0_dp
278 REAL(kind=dp) :: econfinement = 0.0_dp
279 END TYPE atom_energy_type
280
281!> \brief Information on optimization procedure
282! **************************************************************************************************
284 REAL(kind=dp) :: damping = 0.0_dp
285 REAL(kind=dp) :: eps_scf = 0.0_dp
286 REAL(kind=dp) :: eps_diis = 0.0_dp
287 INTEGER :: max_iter = 0
288 INTEGER :: n_diis = 0
290
291!> \brief Provides all information about an atomic kind
292! **************************************************************************************************
294 INTEGER :: z = 0
295 INTEGER :: zcore = 0
296 LOGICAL :: pp_calc = .false.
297! ZMP adding in type some variables
298 LOGICAL :: do_zmp = .false., doread = .false., read_vxc = .false., dm = .false.
299 CHARACTER(LEN=default_string_length) :: ext_file = "", ext_vxc_file = "", &
300 zmp_restart_file = ""
301!
302 INTEGER :: method_type = do_rks_atom
303 INTEGER :: relativistic = do_nonrel_atom
304 INTEGER :: coulomb_integral_type = do_analytic
305 INTEGER :: exchange_integral_type = do_analytic
306! ZMP
307 REAL(kind=dp) :: lambda = 0.0_dp
308 REAL(kind=dp) :: rho_diff_integral = 0.0_dp
309 REAL(kind=dp) :: weight = 0.0_dp, zmpgrid_tol = 0.0_dp, zmpvxcgrid_tol = 0.0_dp
310!
311 TYPE(atom_basis_type), POINTER :: basis => null()
312 TYPE(atom_potential_type), POINTER :: potential => null()
313 TYPE(atom_state), POINTER :: state => null()
314 TYPE(atom_integrals), POINTER :: integrals => null()
315 TYPE(atom_orbitals), POINTER :: orbitals => null()
316 TYPE(atom_energy_type) :: energy = atom_energy_type()
318 TYPE(section_vals_type), POINTER :: xc_section => null(), zmp_section => null()
319 TYPE(opmat_type), POINTER :: fmat => null()
320 TYPE(atom_hfx_type) :: hfx_pot = atom_hfx_type()
321 END TYPE atom_type
322! **************************************************************************************************
324 TYPE(atom_type), POINTER :: atom => null()
325 END TYPE atom_p_type
326
327 PUBLIC :: lmat
331 PUBLIC :: atom_optimization_type
332 PUBLIC :: atom_compare_grids
338 PUBLIC :: clementi_geobas
343 PUBLIC :: setup_hf_section
344
345! **************************************************************************************************
346
347CONTAINS
348
349! **************************************************************************************************
350!> \brief Initialize the basis for the atomic code
351!> \param basis ...
352!> \param basis_section ...
353!> \param zval ...
354!> \param btyp ...
355!> \note Highly accurate relativistic universal Gaussian basis set: Dirac-Fock-Coulomb calculations
356!> for atomic systems up to nobelium
357!> J. Chem. Phys. 101, 6829 (1994); DOI:10.1063/1.468311
358!> G. L. Malli and A. B. F. Da Silva
359!> Department of Chemistry, Simon Fraser University, Burnaby, B.C., Canada
360!> Yasuyuki Ishikawa
361!> Department of Chemistry, University of Puerto Rico, San Juan, Puerto Rico
362!>
363!> A universal Gaussian basis set is developed that leads to relativistic Dirac-Fock SCF energies
364!> of comparable accuracy as that obtained by the accurate numerical finite-difference method
365!> (GRASP2 package) [J. Phys. B 25, 1 (1992)]. The Gaussian-type functions of our universal basis
366!> set satisfy the relativistic boundary conditions associated with the finite nuclear model for a
367!> finite speed of light and conform to the so-called kinetic balance at the nonrelativistic limit.
368!> We attribute the exceptionally high accuracy obtained in our calculations to the fact that the
369!> representation of the relativistic dynamics of an electron in a spherical ball finite nucleus
370!> near the origin in terms of our universal Gaussian basis set is as accurate as that provided by
371!> the numerical finite-difference method. Results of the Dirac-Fock-Coulomb energies for a number
372!> of atoms up to No (Z=102) and some negative ions are presented and compared with the recent
373!> results obtained with the numerical finite-difference method and geometrical Gaussian basis sets
374!> by Parpia, Mohanty, and Clementi [J. Phys. B 25, 1 (1992)]. The accuracy of our calculations is
375!> estimated to be within a few parts in 109 for all the atomic systems studied.
376! **************************************************************************************************
377 SUBROUTINE init_atom_basis(basis, basis_section, zval, btyp)
378 TYPE(atom_basis_type), INTENT(INOUT) :: basis
379 TYPE(section_vals_type), POINTER :: basis_section
380 INTEGER, INTENT(IN) :: zval
381 CHARACTER(LEN=2) :: btyp
382
383 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_atom_basis'
384 INTEGER, PARAMETER :: nua = 40, nup = 16
385 REAL(kind=dp), DIMENSION(nua), PARAMETER :: ugbs = [0.007299_dp, 0.013705_dp, 0.025733_dp, &
386 0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
387 3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
388 174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
389 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
390 94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
391 2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
392 27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
393 341890134.751331_dp]
394
395 CHARACTER(LEN=default_string_length) :: basis_fn, basis_name
396 INTEGER :: basistype, handle, i, j, k, l, ll, m, &
397 ngp, nl, nr, nu, quadtype
398 INTEGER, DIMENSION(0:lmat) :: starti
399 INTEGER, DIMENSION(:), POINTER :: nqm, num_gto, num_slater, sindex
400 REAL(kind=dp) :: al, amax, aval, cval, ear, pf, rk
401 REAL(kind=dp), DIMENSION(:), POINTER :: expo
402 TYPE(section_vals_type), POINTER :: gto_basis_section
403
404 CALL timeset(routinen, handle)
405
406 ! btyp = AE : standard all-electron basis
407 ! btyp = PP : standard pseudopotential basis
408 ! btyp = AA : high accuracy all-electron basis
409 ! btyp = AP : high accuracy pseudopotential basis
410
411 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
412 ! get information on quadrature type and number of grid points
413 ! allocate and initialize the atomic grid
414 CALL allocate_grid_atom(basis%grid)
415 CALL section_vals_val_get(basis_section, "QUADRATURE", i_val=quadtype)
416 CALL section_vals_val_get(basis_section, "GRID_POINTS", i_val=ngp)
417 IF (ngp <= 0) &
418 cpabort("The number of radial grid points must be greater than zero.")
419 CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
420 basis%grid%nr = ngp
421 basis%geometrical = .false.
422 basis%aval = 0._dp
423 basis%cval = 0._dp
424 basis%start = 0
425
426 CALL section_vals_val_get(basis_section, "BASIS_TYPE", i_val=basistype)
427 CALL section_vals_val_get(basis_section, "EPS_EIGENVALUE", r_val=basis%eps_eig)
428 SELECT CASE (basistype)
429 CASE (gaussian)
430 basis%basis_type = gto_basis
431 NULLIFY (num_gto)
432 CALL section_vals_val_get(basis_section, "NUM_GTO", i_vals=num_gto)
433 IF (num_gto(1) < 1) THEN
434 ! use default basis
435 IF (btyp == "AE") THEN
436 nu = nua
437 ELSE IF (btyp == "PP") THEN
438 nu = nup
439 ELSE
440 nu = nua
441 END IF
442 basis%nbas = nu
443 basis%nprim = nu
444 ALLOCATE (basis%am(nu, 0:lmat))
445 DO i = 0, lmat
446 basis%am(1:nu, i) = ugbs(1:nu)
447 END DO
448 ELSE
449 basis%nbas = 0
450 DO i = 1, SIZE(num_gto)
451 basis%nbas(i - 1) = num_gto(i)
452 END DO
453 basis%nprim = basis%nbas
454 m = maxval(basis%nbas)
455 ALLOCATE (basis%am(m, 0:lmat))
456 basis%am = 0._dp
457 DO l = 0, lmat
458 IF (basis%nbas(l) > 0) THEN
459 NULLIFY (expo)
460 SELECT CASE (l)
461 CASE (0)
462 CALL section_vals_val_get(basis_section, "S_EXPONENTS", r_vals=expo)
463 CASE (1)
464 CALL section_vals_val_get(basis_section, "P_EXPONENTS", r_vals=expo)
465 CASE (2)
466 CALL section_vals_val_get(basis_section, "D_EXPONENTS", r_vals=expo)
467 CASE (3)
468 CALL section_vals_val_get(basis_section, "F_EXPONENTS", r_vals=expo)
469 CASE DEFAULT
470 cpabort("Invalid angular quantum number l found for Gaussian basis set")
471 END SELECT
472 cpassert(SIZE(expo) >= basis%nbas(l))
473 DO i = 1, basis%nbas(l)
474 basis%am(i, l) = expo(i)
475 END DO
476 END IF
477 END DO
478 END IF
479 ! initialize basis function on a radial grid
480 nr = basis%grid%nr
481 m = maxval(basis%nbas)
482 ALLOCATE (basis%bf(nr, m, 0:lmat))
483 ALLOCATE (basis%dbf(nr, m, 0:lmat))
484 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
485 basis%bf = 0._dp
486 basis%dbf = 0._dp
487 basis%ddbf = 0._dp
488 DO l = 0, lmat
489 DO i = 1, basis%nbas(l)
490 al = basis%am(i, l)
491 DO k = 1, nr
492 rk = basis%grid%rad(k)
493 ear = exp(-al*basis%grid%rad(k)**2)
494 basis%bf(k, i, l) = rk**l*ear
495 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
496 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
497 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
498 END DO
499 END DO
500 END DO
501 CASE (geometrical_gto)
502 basis%basis_type = gto_basis
503 NULLIFY (num_gto)
504 CALL section_vals_val_get(basis_section, "NUM_GTO", i_vals=num_gto)
505 IF (num_gto(1) < 1) THEN
506 IF (btyp == "AE") THEN
507 ! use the Clementi extra large basis
508 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
509 ELSE IF (btyp == "PP") THEN
510 ! use the Clementi extra large basis
511 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
512 ELSE IF (btyp == "AA") THEN
513 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
514 amax = cval**(basis%nbas(0) - 1)
515 basis%nbas(0) = nint((log(amax)/log(1.6_dp)))
516 cval = 1.6_dp
517 starti = 0
518 basis%nbas(1) = basis%nbas(0) - 4
519 basis%nbas(2) = basis%nbas(0) - 8
520 basis%nbas(3) = basis%nbas(0) - 12
521 IF (lmat > 3) basis%nbas(4:lmat) = 0
522 ELSE IF (btyp == "AP") THEN
523 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
524 amax = 500._dp/aval
525 basis%nbas = nint((log(amax)/log(1.6_dp)))
526 cval = 1.6_dp
527 starti = 0
528 ELSE
529 ! use the Clementi extra large basis
530 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
531 END IF
532 basis%nprim = basis%nbas
533 ELSE
534 basis%nbas = 0
535 DO i = 1, SIZE(num_gto)
536 basis%nbas(i - 1) = num_gto(i)
537 END DO
538 basis%nprim = basis%nbas
539 NULLIFY (sindex)
540 CALL section_vals_val_get(basis_section, "START_INDEX", i_vals=sindex)
541 starti = 0
542 DO i = 1, SIZE(sindex)
543 starti(i - 1) = sindex(i)
544 cpassert(sindex(i) >= 0)
545 END DO
546 CALL section_vals_val_get(basis_section, "GEOMETRICAL_FACTOR", r_val=cval)
547 CALL section_vals_val_get(basis_section, "GEO_START_VALUE", r_val=aval)
548 END IF
549 m = maxval(basis%nbas)
550 ALLOCATE (basis%am(m, 0:lmat))
551 basis%am = 0._dp
552 DO l = 0, lmat
553 DO i = 1, basis%nbas(l)
554 ll = i - 1 + starti(l)
555 basis%am(i, l) = aval*cval**(ll)
556 END DO
557 END DO
558
559 basis%geometrical = .true.
560 basis%aval = aval
561 basis%cval = cval
562 basis%start = starti
563
564 ! initialize basis function on a radial grid
565 nr = basis%grid%nr
566 m = maxval(basis%nbas)
567 ALLOCATE (basis%bf(nr, m, 0:lmat))
568 ALLOCATE (basis%dbf(nr, m, 0:lmat))
569 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
570 basis%bf = 0._dp
571 basis%dbf = 0._dp
572 basis%ddbf = 0._dp
573 DO l = 0, lmat
574 DO i = 1, basis%nbas(l)
575 al = basis%am(i, l)
576 DO k = 1, nr
577 rk = basis%grid%rad(k)
578 ear = exp(-al*basis%grid%rad(k)**2)
579 basis%bf(k, i, l) = rk**l*ear
580 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
581 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
582 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
583 END DO
584 END DO
585 END DO
586 CASE (contracted_gto)
587 basis%basis_type = cgto_basis
588 CALL section_vals_val_get(basis_section, "BASIS_SET_FILE_NAME", c_val=basis_fn)
589 CALL section_vals_val_get(basis_section, "BASIS_SET", c_val=basis_name)
590 gto_basis_section => section_vals_get_subs_vals(basis_section, "BASIS")
591 CALL read_basis_set(ptable(zval)%symbol, basis, basis_name, basis_fn, &
592 gto_basis_section)
593
594 ! initialize basis function on a radial grid
595 nr = basis%grid%nr
596 m = maxval(basis%nbas)
597 ALLOCATE (basis%bf(nr, m, 0:lmat))
598 ALLOCATE (basis%dbf(nr, m, 0:lmat))
599 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
600 basis%bf = 0._dp
601 basis%dbf = 0._dp
602 basis%ddbf = 0._dp
603 DO l = 0, lmat
604 DO i = 1, basis%nprim(l)
605 al = basis%am(i, l)
606 DO k = 1, nr
607 rk = basis%grid%rad(k)
608 ear = exp(-al*basis%grid%rad(k)**2)
609 DO j = 1, basis%nbas(l)
610 basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
611 basis%dbf(k, j, l) = basis%dbf(k, j, l) &
612 + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
613 basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
614 (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))* &
615 ear*basis%cm(i, j, l)
616 END DO
617 END DO
618 END DO
619 END DO
620 CASE (slater)
621 basis%basis_type = sto_basis
622 NULLIFY (num_slater)
623 CALL section_vals_val_get(basis_section, "NUM_SLATER", i_vals=num_slater)
624 IF (num_slater(1) < 1) THEN
625 cpabort("Invalid number (less than 1) Slater-type functions found.")
626 ELSE
627 basis%nbas = 0
628 DO i = 1, SIZE(num_slater)
629 basis%nbas(i - 1) = num_slater(i)
630 END DO
631 basis%nprim = basis%nbas
632 m = maxval(basis%nbas)
633 ALLOCATE (basis%as(m, 0:lmat), basis%ns(m, 0:lmat))
634 basis%as = 0.0_dp
635 basis%ns = 0
636 DO l = 0, lmat
637 IF (basis%nbas(l) > 0) THEN
638 NULLIFY (expo)
639 SELECT CASE (l)
640 CASE (0)
641 CALL section_vals_val_get(basis_section, "S_EXPONENTS", r_vals=expo)
642 CASE (1)
643 CALL section_vals_val_get(basis_section, "P_EXPONENTS", r_vals=expo)
644 CASE (2)
645 CALL section_vals_val_get(basis_section, "D_EXPONENTS", r_vals=expo)
646 CASE (3)
647 CALL section_vals_val_get(basis_section, "F_EXPONENTS", r_vals=expo)
648 CASE DEFAULT
649 cpabort("Invalid angular quantum number l found for Slater basis set")
650 END SELECT
651 cpassert(SIZE(expo) >= basis%nbas(l))
652 DO i = 1, basis%nbas(l)
653 basis%as(i, l) = expo(i)
654 END DO
655 NULLIFY (nqm)
656 SELECT CASE (l)
657 CASE (0)
658 CALL section_vals_val_get(basis_section, "S_QUANTUM_NUMBERS", i_vals=nqm)
659 CASE (1)
660 CALL section_vals_val_get(basis_section, "P_QUANTUM_NUMBERS", i_vals=nqm)
661 CASE (2)
662 CALL section_vals_val_get(basis_section, "D_QUANTUM_NUMBERS", i_vals=nqm)
663 CASE (3)
664 CALL section_vals_val_get(basis_section, "F_QUANTUM_NUMBERS", i_vals=nqm)
665 CASE DEFAULT
666 cpabort("Invalid angular quantum number l found for Slater basis set")
667 END SELECT
668 cpassert(SIZE(nqm) >= basis%nbas(l))
669 DO i = 1, basis%nbas(l)
670 basis%ns(i, l) = nqm(i)
671 END DO
672 END IF
673 END DO
674 END IF
675 ! initialize basis function on a radial grid
676 nr = basis%grid%nr
677 m = maxval(basis%nbas)
678 ALLOCATE (basis%bf(nr, m, 0:lmat))
679 ALLOCATE (basis%dbf(nr, m, 0:lmat))
680 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
681 basis%bf = 0._dp
682 basis%dbf = 0._dp
683 basis%ddbf = 0._dp
684 DO l = 0, lmat
685 DO i = 1, basis%nbas(l)
686 al = basis%as(i, l)
687 nl = basis%ns(i, l)
688 pf = (2._dp*al)**nl*sqrt(2._dp*al/fac(2*nl))
689 DO k = 1, nr
690 rk = basis%grid%rad(k)
691 ear = rk**(nl - 1)*exp(-al*rk)
692 basis%bf(k, i, l) = pf*ear
693 basis%dbf(k, i, l) = pf*(real(nl - 1, dp)/rk - al)*ear
694 basis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1), dp)/rk/rk &
695 - al*real(2*(nl - 1), dp)/rk + al*al)*ear
696 END DO
697 END DO
698 END DO
699 CASE (numerical)
700 basis%basis_type = num_basis
701 cpabort("Numerical basis set type not yet implemented.")
702 CASE DEFAULT
703 cpabort("Unknown basis set type specified. Check the code!")
704 END SELECT
705
706 CALL timestop(handle)
707
708 END SUBROUTINE init_atom_basis
709
710! **************************************************************************************************
711!> \brief ...
712!> \param basis ...
713! **************************************************************************************************
715 TYPE(atom_basis_type), INTENT(INOUT) :: basis
716
717 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_atom_basis_default_pp'
718 INTEGER, PARAMETER :: nua = 40, nup = 20
719 REAL(kind=dp), DIMENSION(nua), PARAMETER :: ugbs = [0.007299_dp, 0.013705_dp, 0.025733_dp, &
720 0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
721 3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
722 174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
723 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
724 94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
725 2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
726 27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
727 341890134.751331_dp]
728
729 INTEGER :: handle, i, k, l, m, ngp, nr, nu, quadtype
730 REAL(kind=dp) :: al, ear, rk
731
732 CALL timeset(routinen, handle)
733
734 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
735
736 ! Allocate and initialize the atomic grid
737 NULLIFY (basis%grid)
738 CALL allocate_grid_atom(basis%grid)
739 quadtype = do_gapw_log
740 ngp = 500
741 CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
742 basis%grid%nr = ngp
743 basis%geometrical = .false.
744 basis%aval = 0._dp
745 basis%cval = 0._dp
746 basis%start = 0
747 basis%eps_eig = 1.e-12_dp
748
749 basis%basis_type = gto_basis
750 nu = nup
751 basis%nbas = nu
752 basis%nprim = nu
753 ALLOCATE (basis%am(nu, 0:lmat))
754 DO i = 0, lmat
755 basis%am(1:nu, i) = ugbs(1:nu)
756 END DO
757 ! initialize basis function on a radial grid
758 nr = basis%grid%nr
759 m = maxval(basis%nbas)
760 ALLOCATE (basis%bf(nr, m, 0:lmat))
761 ALLOCATE (basis%dbf(nr, m, 0:lmat))
762 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
763 basis%bf = 0._dp
764 basis%dbf = 0._dp
765 basis%ddbf = 0._dp
766 DO l = 0, lmat
767 DO i = 1, basis%nbas(l)
768 al = basis%am(i, l)
769 DO k = 1, nr
770 rk = basis%grid%rad(k)
771 ear = exp(-al*basis%grid%rad(k)**2)
772 basis%bf(k, i, l) = rk**l*ear
773 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
774 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
775 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
776 END DO
777 END DO
778 END DO
779
780 CALL timestop(handle)
781
782 END SUBROUTINE init_atom_basis_default_pp
783
784! **************************************************************************************************
785!> \brief ...
786!> \param basis ...
787!> \param gbasis ...
788!> \param r ...
789!> \param rab ...
790! **************************************************************************************************
791 SUBROUTINE atom_basis_gridrep(basis, gbasis, r, rab)
792 TYPE(atom_basis_type), INTENT(IN) :: basis
793 TYPE(atom_basis_type), INTENT(INOUT) :: gbasis
794 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, rab
795
796 INTEGER :: i, j, k, l, m, n1, n2, n3, ngp, nl, nr, &
797 quadtype
798 REAL(kind=dp) :: al, ear, pf, rk
799
800 NULLIFY (gbasis%am, gbasis%cm, gbasis%as, gbasis%ns, gbasis%bf, gbasis%dbf, gbasis%ddbf)
801
802 ! copy basis info
803 gbasis%basis_type = basis%basis_type
804 gbasis%nbas(0:lmat) = basis%nbas(0:lmat)
805 gbasis%nprim(0:lmat) = basis%nprim(0:lmat)
806 IF (ASSOCIATED(basis%am)) THEN
807 n1 = SIZE(basis%am, 1)
808 n2 = SIZE(basis%am, 2)
809 ALLOCATE (gbasis%am(n1, 0:n2 - 1))
810 gbasis%am = basis%am
811 END IF
812 IF (ASSOCIATED(basis%cm)) THEN
813 n1 = SIZE(basis%cm, 1)
814 n2 = SIZE(basis%cm, 2)
815 n3 = SIZE(basis%cm, 3)
816 ALLOCATE (gbasis%cm(n1, n2, 0:n3 - 1))
817 gbasis%cm = basis%cm
818 END IF
819 IF (ASSOCIATED(basis%as)) THEN
820 n1 = SIZE(basis%as, 1)
821 n2 = SIZE(basis%as, 2)
822 ALLOCATE (gbasis%as(n1, 0:n2 - 1))
823 gbasis%as = basis%as
824 END IF
825 IF (ASSOCIATED(basis%ns)) THEN
826 n1 = SIZE(basis%ns, 1)
827 n2 = SIZE(basis%ns, 2)
828 ALLOCATE (gbasis%ns(n1, 0:n2 - 1))
829 gbasis%ns = basis%ns
830 END IF
831 gbasis%eps_eig = basis%eps_eig
832 gbasis%geometrical = basis%geometrical
833 gbasis%aval = basis%aval
834 gbasis%cval = basis%cval
835 gbasis%start(0:lmat) = basis%start(0:lmat)
836
837 ! get information on quadrature type and number of grid points
838 ! allocate and initialize the atomic grid
839 NULLIFY (gbasis%grid)
840 CALL allocate_grid_atom(gbasis%grid)
841 ngp = SIZE(r)
842 quadtype = do_gapw_log
843 IF (ngp <= 0) &
844 cpabort("The number of radial grid points must be greater than zero.")
845 CALL create_grid_atom(gbasis%grid, ngp, 1, 1, 0, quadtype)
846 gbasis%grid%nr = ngp
847 gbasis%grid%rad(:) = r(:)
848 gbasis%grid%rad2(:) = r(:)*r(:)
849 gbasis%grid%wr(:) = rab(:)*gbasis%grid%rad2(:)
850
851 ! initialize basis function on a radial grid
852 nr = gbasis%grid%nr
853 m = maxval(gbasis%nbas)
854 ALLOCATE (gbasis%bf(nr, m, 0:lmat))
855 ALLOCATE (gbasis%dbf(nr, m, 0:lmat))
856 ALLOCATE (gbasis%ddbf(nr, m, 0:lmat))
857 gbasis%bf = 0._dp
858 gbasis%dbf = 0._dp
859 gbasis%ddbf = 0._dp
860
861 SELECT CASE (gbasis%basis_type)
862 CASE (gto_basis)
863 DO l = 0, lmat
864 DO i = 1, gbasis%nbas(l)
865 al = gbasis%am(i, l)
866 DO k = 1, nr
867 rk = gbasis%grid%rad(k)
868 ear = exp(-al*gbasis%grid%rad(k)**2)
869 gbasis%bf(k, i, l) = rk**l*ear
870 gbasis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
871 gbasis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
872 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
873 END DO
874 END DO
875 END DO
876 CASE (cgto_basis)
877 DO l = 0, lmat
878 DO i = 1, gbasis%nprim(l)
879 al = gbasis%am(i, l)
880 DO k = 1, nr
881 rk = gbasis%grid%rad(k)
882 ear = exp(-al*gbasis%grid%rad(k)**2)
883 DO j = 1, gbasis%nbas(l)
884 gbasis%bf(k, j, l) = gbasis%bf(k, j, l) + rk**l*ear*gbasis%cm(i, j, l)
885 gbasis%dbf(k, j, l) = gbasis%dbf(k, j, l) &
886 + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*gbasis%cm(i, j, l)
887 gbasis%ddbf(k, j, l) = gbasis%ddbf(k, j, l) + &
888 (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))* &
889 ear*gbasis%cm(i, j, l)
890 END DO
891 END DO
892 END DO
893 END DO
894 CASE (sto_basis)
895 DO l = 0, lmat
896 DO i = 1, gbasis%nbas(l)
897 al = gbasis%as(i, l)
898 nl = gbasis%ns(i, l)
899 pf = (2._dp*al)**nl*sqrt(2._dp*al/fac(2*nl))
900 DO k = 1, nr
901 rk = gbasis%grid%rad(k)
902 ear = rk**(nl - 1)*exp(-al*rk)
903 gbasis%bf(k, i, l) = pf*ear
904 gbasis%dbf(k, i, l) = pf*(real(nl - 1, dp)/rk - al)*ear
905 gbasis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1), dp)/rk/rk &
906 - al*real(2*(nl - 1), dp)/rk + al*al)*ear
907 END DO
908 END DO
909 END DO
910 CASE (num_basis)
911 gbasis%basis_type = num_basis
912 cpabort("Numerical basis set type not yet implemented.")
913 CASE DEFAULT
914 cpabort("Unknown basis set type specified. Check the code!")
915 END SELECT
916
917 END SUBROUTINE atom_basis_gridrep
918
919! **************************************************************************************************
920!> \brief ...
921!> \param basis ...
922! **************************************************************************************************
923 SUBROUTINE release_atom_basis(basis)
924 TYPE(atom_basis_type), INTENT(INOUT) :: basis
925
926 IF (ASSOCIATED(basis%am)) THEN
927 DEALLOCATE (basis%am)
928 END IF
929 IF (ASSOCIATED(basis%cm)) THEN
930 DEALLOCATE (basis%cm)
931 END IF
932 IF (ASSOCIATED(basis%as)) THEN
933 DEALLOCATE (basis%as)
934 END IF
935 IF (ASSOCIATED(basis%ns)) THEN
936 DEALLOCATE (basis%ns)
937 END IF
938 IF (ASSOCIATED(basis%bf)) THEN
939 DEALLOCATE (basis%bf)
940 END IF
941 IF (ASSOCIATED(basis%dbf)) THEN
942 DEALLOCATE (basis%dbf)
943 END IF
944 IF (ASSOCIATED(basis%ddbf)) THEN
945 DEALLOCATE (basis%ddbf)
946 END IF
947
948 CALL deallocate_grid_atom(basis%grid)
949
950 END SUBROUTINE release_atom_basis
951! **************************************************************************************************
952
953! **************************************************************************************************
954!> \brief ...
955!> \param atom ...
956! **************************************************************************************************
957 SUBROUTINE create_atom_type(atom)
958 TYPE(atom_type), POINTER :: atom
959
960 cpassert(.NOT. ASSOCIATED(atom))
961
962 ALLOCATE (atom)
963
964 NULLIFY (atom%zmp_section)
965 NULLIFY (atom%xc_section)
966 NULLIFY (atom%fmat)
967 atom%do_zmp = .false.
968 atom%doread = .false.
969 atom%read_vxc = .false.
970 atom%dm = .false.
971 atom%hfx_pot%scale_coulomb = 0.0_dp
972 atom%hfx_pot%scale_longrange = 0.0_dp
973 atom%hfx_pot%omega = 0.0_dp
974
975 END SUBROUTINE create_atom_type
976
977! **************************************************************************************************
978!> \brief ...
979!> \param atom ...
980! **************************************************************************************************
981 SUBROUTINE release_atom_type(atom)
982 TYPE(atom_type), POINTER :: atom
983
984 cpassert(ASSOCIATED(atom))
985
986 NULLIFY (atom%basis)
987 NULLIFY (atom%integrals)
988 IF (ASSOCIATED(atom%state)) THEN
989 DEALLOCATE (atom%state)
990 END IF
991 IF (ASSOCIATED(atom%orbitals)) THEN
992 CALL release_atom_orbs(atom%orbitals)
993 END IF
994
995 IF (ASSOCIATED(atom%fmat)) CALL release_opmat(atom%fmat)
996
997 DEALLOCATE (atom)
998
999 END SUBROUTINE release_atom_type
1000
1001! ZMP adding input variables in subroutine do_zmp,doread,read_vxc,method_type
1002! **************************************************************************************************
1003!> \brief ...
1004!> \param atom ...
1005!> \param basis ...
1006!> \param state ...
1007!> \param integrals ...
1008!> \param orbitals ...
1009!> \param potential ...
1010!> \param zcore ...
1011!> \param pp_calc ...
1012!> \param do_zmp ...
1013!> \param doread ...
1014!> \param read_vxc ...
1015!> \param method_type ...
1016!> \param relativistic ...
1017!> \param coulomb_integral_type ...
1018!> \param exchange_integral_type ...
1019!> \param fmat ...
1020! **************************************************************************************************
1021 SUBROUTINE set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, &
1022 read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
1023 TYPE(atom_type), POINTER :: atom
1024 TYPE(atom_basis_type), OPTIONAL, POINTER :: basis
1025 TYPE(atom_state), OPTIONAL, POINTER :: state
1026 TYPE(atom_integrals), OPTIONAL, POINTER :: integrals
1027 TYPE(atom_orbitals), OPTIONAL, POINTER :: orbitals
1028 TYPE(atom_potential_type), OPTIONAL, POINTER :: potential
1029 INTEGER, INTENT(IN), OPTIONAL :: zcore
1030 LOGICAL, INTENT(IN), OPTIONAL :: pp_calc, do_zmp, doread, read_vxc
1031 INTEGER, INTENT(IN), OPTIONAL :: method_type, relativistic, &
1032 coulomb_integral_type, &
1033 exchange_integral_type
1034 TYPE(opmat_type), OPTIONAL, POINTER :: fmat
1035
1036 cpassert(ASSOCIATED(atom))
1037
1038 IF (PRESENT(basis)) atom%basis => basis
1039 IF (PRESENT(state)) atom%state => state
1040 IF (PRESENT(integrals)) atom%integrals => integrals
1041 IF (PRESENT(orbitals)) atom%orbitals => orbitals
1042 IF (PRESENT(potential)) atom%potential => potential
1043 IF (PRESENT(zcore)) atom%zcore = zcore
1044 IF (PRESENT(pp_calc)) atom%pp_calc = pp_calc
1045! ZMP assigning variable values if present
1046 IF (PRESENT(do_zmp)) atom%do_zmp = do_zmp
1047 IF (PRESENT(doread)) atom%doread = doread
1048 IF (PRESENT(read_vxc)) atom%read_vxc = read_vxc
1049
1050 IF (PRESENT(method_type)) atom%method_type = method_type
1051 IF (PRESENT(relativistic)) atom%relativistic = relativistic
1052 IF (PRESENT(coulomb_integral_type)) atom%coulomb_integral_type = coulomb_integral_type
1053 IF (PRESENT(exchange_integral_type)) atom%exchange_integral_type = exchange_integral_type
1054
1055 IF (PRESENT(fmat)) THEN
1056 IF (ASSOCIATED(atom%fmat)) CALL release_opmat(atom%fmat)
1057 atom%fmat => fmat
1058 END IF
1059
1060 END SUBROUTINE set_atom
1061
1062! **************************************************************************************************
1063!> \brief ...
1064!> \param orbs ...
1065!> \param mbas ...
1066!> \param mo ...
1067! **************************************************************************************************
1068 SUBROUTINE create_atom_orbs(orbs, mbas, mo)
1069 TYPE(atom_orbitals), POINTER :: orbs
1070 INTEGER, INTENT(IN) :: mbas, mo
1071
1072 cpassert(.NOT. ASSOCIATED(orbs))
1073
1074 ALLOCATE (orbs)
1075
1076 ALLOCATE (orbs%wfn(mbas, mo, 0:lmat), orbs%wfna(mbas, mo, 0:lmat), orbs%wfnb(mbas, mo, 0:lmat))
1077 orbs%wfn = 0._dp
1078 orbs%wfna = 0._dp
1079 orbs%wfnb = 0._dp
1080
1081 ALLOCATE (orbs%pmat(mbas, mbas, 0:lmat), orbs%pmata(mbas, mbas, 0:lmat), orbs%pmatb(mbas, mbas, 0:lmat))
1082 orbs%pmat = 0._dp
1083 orbs%pmata = 0._dp
1084 orbs%pmatb = 0._dp
1085
1086 ALLOCATE (orbs%ener(mo, 0:lmat), orbs%enera(mo, 0:lmat), orbs%enerb(mo, 0:lmat))
1087 orbs%ener = 0._dp
1088 orbs%enera = 0._dp
1089 orbs%enerb = 0._dp
1090
1091 ALLOCATE (orbs%refene(mo, 0:lmat, 2), orbs%refchg(mo, 0:lmat, 2), orbs%refnod(mo, 0:lmat, 2))
1092 orbs%refene = 0._dp
1093 orbs%refchg = 0._dp
1094 orbs%refnod = 0._dp
1095 ALLOCATE (orbs%wrefene(mo, 0:lmat, 2), orbs%wrefchg(mo, 0:lmat, 2), orbs%wrefnod(mo, 0:lmat, 2))
1096 orbs%wrefene = 0._dp
1097 orbs%wrefchg = 0._dp
1098 orbs%wrefnod = 0._dp
1099 ALLOCATE (orbs%crefene(mo, 0:lmat, 2), orbs%crefchg(mo, 0:lmat, 2), orbs%crefnod(mo, 0:lmat, 2))
1100 orbs%crefene = 0._dp
1101 orbs%crefchg = 0._dp
1102 orbs%crefnod = 0._dp
1103 ALLOCATE (orbs%rcmax(mo, 0:lmat, 2))
1104 orbs%rcmax = 0._dp
1105 ALLOCATE (orbs%wpsir0(mo, 2), orbs%tpsir0(mo, 2))
1106 orbs%wpsir0 = 0._dp
1107 orbs%tpsir0 = 0._dp
1108 ALLOCATE (orbs%reftype(mo, 0:lmat, 2))
1109 orbs%reftype = "XX"
1110
1111 END SUBROUTINE create_atom_orbs
1112
1113! **************************************************************************************************
1114!> \brief ...
1115!> \param orbs ...
1116! **************************************************************************************************
1117 SUBROUTINE release_atom_orbs(orbs)
1118 TYPE(atom_orbitals), POINTER :: orbs
1119
1120 cpassert(ASSOCIATED(orbs))
1121
1122 IF (ASSOCIATED(orbs%wfn)) THEN
1123 DEALLOCATE (orbs%wfn, orbs%wfna, orbs%wfnb)
1124 END IF
1125 IF (ASSOCIATED(orbs%pmat)) THEN
1126 DEALLOCATE (orbs%pmat, orbs%pmata, orbs%pmatb)
1127 END IF
1128 IF (ASSOCIATED(orbs%ener)) THEN
1129 DEALLOCATE (orbs%ener, orbs%enera, orbs%enerb)
1130 END IF
1131 IF (ASSOCIATED(orbs%refene)) THEN
1132 DEALLOCATE (orbs%refene)
1133 END IF
1134 IF (ASSOCIATED(orbs%refchg)) THEN
1135 DEALLOCATE (orbs%refchg)
1136 END IF
1137 IF (ASSOCIATED(orbs%refnod)) THEN
1138 DEALLOCATE (orbs%refnod)
1139 END IF
1140 IF (ASSOCIATED(orbs%wrefene)) THEN
1141 DEALLOCATE (orbs%wrefene)
1142 END IF
1143 IF (ASSOCIATED(orbs%wrefchg)) THEN
1144 DEALLOCATE (orbs%wrefchg)
1145 END IF
1146 IF (ASSOCIATED(orbs%wrefnod)) THEN
1147 DEALLOCATE (orbs%wrefnod)
1148 END IF
1149 IF (ASSOCIATED(orbs%crefene)) THEN
1150 DEALLOCATE (orbs%crefene)
1151 END IF
1152 IF (ASSOCIATED(orbs%crefchg)) THEN
1153 DEALLOCATE (orbs%crefchg)
1154 END IF
1155 IF (ASSOCIATED(orbs%crefnod)) THEN
1156 DEALLOCATE (orbs%crefnod)
1157 END IF
1158 IF (ASSOCIATED(orbs%rcmax)) THEN
1159 DEALLOCATE (orbs%rcmax)
1160 END IF
1161 IF (ASSOCIATED(orbs%wpsir0)) THEN
1162 DEALLOCATE (orbs%wpsir0)
1163 END IF
1164 IF (ASSOCIATED(orbs%tpsir0)) THEN
1165 DEALLOCATE (orbs%tpsir0)
1166 END IF
1167 IF (ASSOCIATED(orbs%reftype)) THEN
1168 DEALLOCATE (orbs%reftype)
1169 END IF
1170
1171 DEALLOCATE (orbs)
1172
1173 END SUBROUTINE release_atom_orbs
1174
1175! **************************************************************************************************
1176!> \brief ...
1177!> \param hf_frac ...
1178!> \param do_hfx ...
1179!> \param atom ...
1180!> \param xc_section ...
1181!> \param extype ...
1182! **************************************************************************************************
1183 SUBROUTINE setup_hf_section(hf_frac, do_hfx, atom, xc_section, extype)
1184 REAL(kind=dp), INTENT(OUT) :: hf_frac
1185 LOGICAL, INTENT(OUT) :: do_hfx
1186 TYPE(atom_type), INTENT(IN), POINTER :: atom
1187 TYPE(section_vals_type), POINTER :: xc_section
1188 INTEGER, INTENT(IN) :: extype
1189
1190 INTEGER :: i, j, nr, nu, pot_type
1191 REAL(kind=dp) :: scale_coulomb, scale_longrange
1192 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: abscissa, weights
1193 TYPE(section_vals_type), POINTER :: hf_sub_section, hfx_sections
1194
1195 hf_frac = 0._dp
1196 IF (ASSOCIATED(atom%xc_section)) THEN
1197 xc_section => atom%xc_section
1198 hfx_sections => section_vals_get_subs_vals(xc_section, "HF")
1199 CALL section_vals_get(hfx_sections, explicit=do_hfx)
1200
1201 ! If nothing has been set explicitly, assume a Coulomb potential
1202 atom%hfx_pot%scale_longrange = 0.0_dp
1203 atom%hfx_pot%scale_coulomb = 1.0_dp
1204
1205 IF (do_hfx) THEN
1206 CALL section_vals_val_get(hfx_sections, "FRACTION", r_val=hf_frac)
1207
1208 ! Get potential info
1209 hf_sub_section => section_vals_get_subs_vals(hfx_sections, "INTERACTION_POTENTIAL", i_rep_section=1)
1210 CALL section_vals_val_get(hf_sub_section, "POTENTIAL_TYPE", i_val=pot_type)
1211 CALL section_vals_val_get(hf_sub_section, "OMEGA", r_val=atom%hfx_pot%omega)
1212 CALL section_vals_val_get(hf_sub_section, "SCALE_COULOMB", r_val=scale_coulomb)
1213 CALL section_vals_val_get(hf_sub_section, "SCALE_LONGRANGE", r_val=scale_longrange)
1214
1215 ! Setup atomic hfx potential
1216 SELECT CASE (pot_type)
1217 CASE DEFAULT
1218 cpwarn("Potential not implemented, use Coulomb instead!")
1220 atom%hfx_pot%scale_longrange = 0.0_dp
1221 atom%hfx_pot%scale_coulomb = scale_coulomb
1222 CASE (do_potential_long)
1223 atom%hfx_pot%scale_coulomb = 0.0_dp
1224 atom%hfx_pot%scale_longrange = scale_longrange
1225 CASE (do_potential_short)
1226 atom%hfx_pot%scale_coulomb = 1.0_dp
1227 atom%hfx_pot%scale_longrange = -1.0_dp
1228 CASE (do_potential_mix_cl)
1229 atom%hfx_pot%scale_coulomb = scale_coulomb
1230 atom%hfx_pot%scale_longrange = scale_longrange
1231 END SELECT
1232 END IF
1233
1234 ! Check whether extype is supported
1235 IF (atom%hfx_pot%scale_longrange /= 0.0_dp .AND. extype /= do_numeric .AND. extype /= do_semi_analytic) THEN
1236 cpabort("Only numerical and semi-analytic lrHF exchange available!")
1237 END IF
1238
1239 IF (atom%hfx_pot%scale_longrange /= 0.0_dp .AND. extype == do_numeric .AND. .NOT. ALLOCATED(atom%hfx_pot%kernel)) THEN
1240 CALL cite_reference(limpanuparb2011)
1241
1242 IF (atom%hfx_pot%do_gh) THEN
1243 ! Setup kernel for Ewald operator
1244 ! Because of the high computational costs of its calculation, we precalculate it here
1245 ! Use Gauss-Hermite grid instead of the external grid
1246 ALLOCATE (weights(atom%hfx_pot%nr_gh), abscissa(atom%hfx_pot%nr_gh))
1247 CALL get_gauss_hermite_weights(abscissa, weights, atom%hfx_pot%nr_gh)
1248
1249 nr = atom%basis%grid%nr
1250 ALLOCATE (atom%hfx_pot%kernel(nr, atom%hfx_pot%nr_gh, 0:atom%state%maxl_calc + atom%state%maxl_occ))
1251 atom%hfx_pot%kernel = 0.0_dp
1252 DO nu = 0, atom%state%maxl_calc + atom%state%maxl_occ
1253 DO i = 1, atom%hfx_pot%nr_gh
1254 DO j = 1, nr
1255 atom%hfx_pot%kernel(j, i, nu) = bessel0(2.0_dp*atom%hfx_pot%omega &
1256 *abscissa(i)*atom%basis%grid%rad(j), nu)*sqrt(weights(i))
1257 END DO
1258 END DO
1259 END DO
1260 ELSE
1261 ! Setup kernel for Ewald operator
1262 ! Because of the high computational costs of its calculation, we precalculate it here
1263 ! Choose it symmetric to further reduce the costs
1264 nr = atom%basis%grid%nr
1265 ALLOCATE (atom%hfx_pot%kernel(nr, nr, 0:atom%state%maxl_calc + atom%state%maxl_occ))
1266 atom%hfx_pot%kernel = 0.0_dp
1267 DO nu = 0, atom%state%maxl_calc + atom%state%maxl_occ
1268 DO i = 1, nr
1269 DO j = 1, i
1270 atom%hfx_pot%kernel(j, i, nu) = bessel0(2.0_dp*atom%hfx_pot%omega &
1271 *atom%basis%grid%rad(i)*atom%basis%grid%rad(j), nu)
1272 END DO
1273 END DO
1274 END DO
1275 END IF
1276 END IF
1277 ELSE
1278 NULLIFY (xc_section)
1279 do_hfx = .false.
1280 END IF
1281
1282 END SUBROUTINE setup_hf_section
1283
1284! **************************************************************************************************
1285!> \brief ...
1286!> \param abscissa ...
1287!> \param weights ...
1288!> \param nn ...
1289! **************************************************************************************************
1290 SUBROUTINE get_gauss_hermite_weights(abscissa, weights, nn)
1291 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: abscissa, weights
1292 INTEGER, INTENT(IN) :: nn
1293
1294 INTEGER :: counter, ii, info, liwork, lwork
1295 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1296 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag, subdiag, work
1297 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvec
1298
1299 ! Setup matrix for Golub-Welsch-algorithm to determine roots and weights of Gauss-Hermite quadrature
1300 ! If necessary, one can setup matrices differently for other quadratures
1301 ALLOCATE (work(1), iwork(1), diag(2*nn), subdiag(2*nn - 1), eigenvec(2*nn, 2*nn))
1302 lwork = -1
1303 liwork = -1
1304 diag = 0.0_dp
1305 DO ii = 1, 2*nn - 1
1306 subdiag(ii) = sqrt(real(ii, kind=dp)/2.0_dp)
1307 END DO
1308
1309 ! Get correct size for working matrices
1310 CALL dstevd('V', 2*nn, diag, subdiag, eigenvec, 2*nn, work, lwork, iwork, liwork, info)
1311 IF (info /= 0) THEN
1312 ! This should not happen!
1313 cpabort('Finding size of working matrices failed!')
1314 END IF
1315
1316 ! Setup working matrices with their respective optimal sizes
1317 lwork = int(work(1))
1318 liwork = iwork(1)
1319 DEALLOCATE (work, iwork)
1320 ALLOCATE (work(lwork), iwork(liwork))
1321
1322 ! Perform the actual eigenvalue decomposition
1323 CALL dstevd('V', 2*nn, diag, subdiag, eigenvec, 2*nn, work, lwork, iwork, liwork, info)
1324 IF (info /= 0) THEN
1325 ! This should not happen for the usual values of nn! (Checked for nn = 2000)
1326 cpabort('Eigenvalue decomposition failed!')
1327 END IF
1328
1329 DEALLOCATE (work, iwork, subdiag)
1330
1331 ! Identify positive roots of hermite polynomials (zeros of Hermite polynomials are symmetric wrt the origin)
1332 ! We will only keep the positive roots
1333 counter = 0
1334 DO ii = 1, 2*nn
1335 IF (diag(ii) > 0.0_dp) THEN
1336 counter = counter + 1
1337 abscissa(counter) = diag(ii)
1338 weights(counter) = rootpi*eigenvec(1, ii)**2
1339 END IF
1340 END DO
1341 IF (counter /= nn) THEN
1342 cpabort('Have not found enough or too many zeros!')
1343 END IF
1344
1345 END SUBROUTINE get_gauss_hermite_weights
1346
1347! **************************************************************************************************
1348!> \brief ...
1349!> \param opmat ...
1350!> \param n ...
1351!> \param lmax ...
1352! **************************************************************************************************
1353 SUBROUTINE create_opmat(opmat, n, lmax)
1354 TYPE(opmat_type), POINTER :: opmat
1355 INTEGER, DIMENSION(0:lmat), INTENT(IN) :: n
1356 INTEGER, INTENT(IN), OPTIONAL :: lmax
1357
1358 INTEGER :: lm, m
1359
1360 m = maxval(n)
1361 IF (PRESENT(lmax)) THEN
1362 lm = lmax
1363 ELSE
1364 lm = lmat
1365 END IF
1366
1367 cpassert(.NOT. ASSOCIATED(opmat))
1368
1369 ALLOCATE (opmat)
1370
1371 opmat%n = n
1372 ALLOCATE (opmat%op(m, m, 0:lm))
1373 opmat%op = 0._dp
1374
1375 END SUBROUTINE create_opmat
1376
1377! **************************************************************************************************
1378!> \brief ...
1379!> \param opmat ...
1380! **************************************************************************************************
1381 SUBROUTINE release_opmat(opmat)
1382 TYPE(opmat_type), POINTER :: opmat
1383
1384 cpassert(ASSOCIATED(opmat))
1385
1386 opmat%n = 0
1387 DEALLOCATE (opmat%op)
1388
1389 DEALLOCATE (opmat)
1390
1391 END SUBROUTINE release_opmat
1392
1393! **************************************************************************************************
1394!> \brief ...
1395!> \param opgrid ...
1396!> \param grid ...
1397! **************************************************************************************************
1398 SUBROUTINE create_opgrid(opgrid, grid)
1399 TYPE(opgrid_type), POINTER :: opgrid
1400 TYPE(grid_atom_type), POINTER :: grid
1401
1402 INTEGER :: nr
1403
1404 cpassert(.NOT. ASSOCIATED(opgrid))
1405
1406 ALLOCATE (opgrid)
1407
1408 opgrid%grid => grid
1409
1410 nr = grid%nr
1411
1412 ALLOCATE (opgrid%op(nr))
1413 opgrid%op = 0._dp
1414
1415 END SUBROUTINE create_opgrid
1416
1417! **************************************************************************************************
1418!> \brief ...
1419!> \param opgrid ...
1420! **************************************************************************************************
1421 SUBROUTINE release_opgrid(opgrid)
1422 TYPE(opgrid_type), POINTER :: opgrid
1423
1424 cpassert(ASSOCIATED(opgrid))
1425
1426 NULLIFY (opgrid%grid)
1427 DEALLOCATE (opgrid%op)
1428
1429 DEALLOCATE (opgrid)
1430
1431 END SUBROUTINE release_opgrid
1432
1433! **************************************************************************************************
1434!> \brief ...
1435!> \param zval ...
1436!> \param cval ...
1437!> \param aval ...
1438!> \param ngto ...
1439!> \param ival ...
1440! **************************************************************************************************
1441 SUBROUTINE clementi_geobas(zval, cval, aval, ngto, ival)
1442
1443 INTEGER, INTENT(IN) :: zval
1444 REAL(dp), INTENT(OUT) :: cval, aval
1445 INTEGER, DIMENSION(0:lmat), INTENT(OUT) :: ngto, ival
1446
1447 ngto = 0
1448 ival = 0
1449 cval = 0._dp
1450 aval = 0._dp
1451
1452 SELECT CASE (zval)
1453 CASE (1) ! this is from the general geometrical basis and extended
1454 cval = 2.0_dp
1455 aval = 0.016_dp
1456 ngto(0) = 20
1457 CASE (2)
1458 cval = 2.14774520_dp
1459 aval = 0.04850670_dp
1460 ngto(0) = 20
1461 CASE (3)
1462 cval = 2.08932430_dp
1463 aval = 0.02031060_dp
1464 ngto(0) = 23
1465 CASE (4)
1466 cval = 2.09753060_dp
1467 aval = 0.03207070_dp
1468 ngto(0) = 23
1469 CASE (5)
1470 cval = 2.10343410_dp
1471 aval = 0.03591970_dp
1472 ngto(0) = 23
1473 ngto(1) = 16
1474 CASE (6)
1475 cval = 2.10662820_dp
1476 aval = 0.05292410_dp
1477 ngto(0) = 23
1478 ngto(1) = 16
1479 CASE (7)
1480 cval = 2.13743840_dp
1481 aval = 0.06291970_dp
1482 ngto(0) = 23
1483 ngto(1) = 16
1484 CASE (8)
1485 cval = 2.08687310_dp
1486 aval = 0.08350860_dp
1487 ngto(0) = 23
1488 ngto(1) = 16
1489 CASE (9)
1490 cval = 2.12318180_dp
1491 aval = 0.09899170_dp
1492 ngto(0) = 23
1493 ngto(1) = 16
1494 CASE (10)
1495 cval = 2.13164810_dp
1496 aval = 0.11485350_dp
1497 ngto(0) = 23
1498 ngto(1) = 16
1499 CASE (11)
1500 cval = 2.11413310_dp
1501 aval = 0.00922630_dp
1502 ngto(0) = 26
1503 ngto(1) = 16
1504 ival(1) = 4
1505 CASE (12)
1506 cval = 2.12183620_dp
1507 aval = 0.01215850_dp
1508 ngto(0) = 26
1509 ngto(1) = 16
1510 ival(1) = 4
1511 CASE (13)
1512 cval = 2.06073230_dp
1513 aval = 0.01449350_dp
1514 ngto(0) = 26
1515 ngto(1) = 20
1516 ival(0) = 1
1517 CASE (14)
1518 cval = 2.08563660_dp
1519 aval = 0.01861460_dp
1520 ngto(0) = 26
1521 ngto(1) = 20
1522 ival(0) = 1
1523 CASE (15)
1524 cval = 2.04879270_dp
1525 aval = 0.02147790_dp
1526 ngto(0) = 26
1527 ngto(1) = 20
1528 ival(0) = 1
1529 CASE (16)
1530 cval = 2.06216660_dp
1531 aval = 0.01978920_dp
1532 ngto(0) = 26
1533 ngto(1) = 20
1534 ival(0) = 1
1535 CASE (17)
1536 cval = 2.04628670_dp
1537 aval = 0.02451470_dp
1538 ngto(0) = 26
1539 ngto(1) = 20
1540 ival(0) = 1
1541 CASE (18)
1542 cval = 2.08675200_dp
1543 aval = 0.02635040_dp
1544 ngto(0) = 26
1545 ngto(1) = 20
1546 ival(0) = 1
1547 CASE (19)
1548 cval = 2.02715220_dp
1549 aval = 0.01822040_dp
1550 ngto(0) = 29
1551 ngto(1) = 20
1552 ival(1) = 2
1553 CASE (20)
1554 cval = 2.01465650_dp
1555 aval = 0.01646570_dp
1556 ngto(0) = 29
1557 ngto(1) = 20
1558 ival(1) = 2
1559 CASE (21)
1560 cval = 2.01605240_dp
1561 aval = 0.01254190_dp
1562 ngto(0) = 30
1563 ngto(1) = 20
1564 ngto(2) = 18
1565 ival(1) = 2
1566 CASE (22)
1567 cval = 2.01800000_dp
1568 aval = 0.01195490_dp
1569 ngto(0) = 30
1570 ngto(1) = 21
1571 ngto(2) = 17
1572 ival(1) = 2
1573 ival(2) = 1
1574 CASE (23)
1575 cval = 1.98803560_dp
1576 aval = 0.02492140_dp
1577 ngto(0) = 30
1578 ngto(1) = 21
1579 ngto(2) = 17
1580 ival(1) = 2
1581 ival(2) = 1
1582 CASE (24)
1583 cval = 1.98984000_dp
1584 aval = 0.02568400_dp
1585 ngto(0) = 30
1586 ngto(1) = 21
1587 ngto(2) = 17
1588 ival(1) = 2
1589 ival(2) = 1
1590 CASE (25)
1591 cval = 2.01694380_dp
1592 aval = 0.02664480_dp
1593 ngto(0) = 30
1594 ngto(1) = 21
1595 ngto(2) = 17
1596 ival(1) = 2
1597 ival(2) = 1
1598 CASE (26)
1599 cval = 2.01824090_dp
1600 aval = 0.01355000_dp
1601 ngto(0) = 30
1602 ngto(1) = 21
1603 ngto(2) = 17
1604 ival(1) = 2
1605 ival(2) = 1
1606 CASE (27)
1607 cval = 1.98359400_dp
1608 aval = 0.01702210_dp
1609 ngto(0) = 30
1610 ngto(1) = 21
1611 ngto(2) = 17
1612 ival(1) = 2
1613 ival(2) = 2
1614 CASE (28)
1615 cval = 1.96797340_dp
1616 aval = 0.02163180_dp
1617 ngto(0) = 30
1618 ngto(1) = 22
1619 ngto(2) = 17
1620 ival(1) = 3
1621 ival(2) = 2
1622 CASE (29)
1623 cval = 1.98955180_dp
1624 aval = 0.02304480_dp
1625 ngto(0) = 30
1626 ngto(1) = 20
1627 ngto(2) = 17
1628 ival(1) = 3
1629 ival(2) = 2
1630 CASE (30)
1631 cval = 1.98074320_dp
1632 aval = 0.02754320_dp
1633 ngto(0) = 30
1634 ngto(1) = 21
1635 ngto(2) = 17
1636 ival(1) = 3
1637 ival(2) = 2
1638 CASE (31)
1639 cval = 2.00551070_dp
1640 aval = 0.02005530_dp
1641 ngto(0) = 30
1642 ngto(1) = 23
1643 ngto(2) = 17
1644 ival(0) = 1
1645 ival(2) = 2
1646 CASE (32)
1647 cval = 2.00000030_dp
1648 aval = 0.02003000_dp
1649 ngto(0) = 30
1650 ngto(1) = 24
1651 ngto(2) = 17
1652 ival(0) = 1
1653 ival(2) = 2
1654 CASE (33)
1655 cval = 2.00609100_dp
1656 aval = 0.02055620_dp
1657 ngto(0) = 30
1658 ngto(1) = 23
1659 ngto(2) = 17
1660 ival(0) = 1
1661 ival(2) = 2
1662 CASE (34)
1663 cval = 2.00701000_dp
1664 aval = 0.02230400_dp
1665 ngto(0) = 30
1666 ngto(1) = 24
1667 ngto(2) = 17
1668 ival(0) = 1
1669 ival(2) = 2
1670 CASE (35)
1671 cval = 2.01508710_dp
1672 aval = 0.02685790_dp
1673 ngto(0) = 30
1674 ngto(1) = 24
1675 ngto(2) = 17
1676 ival(0) = 1
1677 ival(2) = 2
1678 CASE (36)
1679 cval = 2.01960430_dp
1680 aval = 0.02960430_dp
1681 ngto(0) = 30
1682 ngto(1) = 24
1683 ngto(2) = 17
1684 ival(0) = 1
1685 ival(2) = 2
1686 CASE (37)
1687 cval = 2.00031000_dp
1688 aval = 0.00768400_dp
1689 ngto(0) = 32
1690 ngto(1) = 25
1691 ngto(2) = 17
1692 ival(0) = 1
1693 ival(1) = 1
1694 ival(2) = 4
1695 CASE (38)
1696 cval = 1.99563960_dp
1697 aval = 0.01401940_dp
1698 ngto(0) = 33
1699 ngto(1) = 24
1700 ngto(2) = 17
1701 ival(1) = 1
1702 ival(2) = 4
1703 CASE (39)
1704 cval = 1.98971210_dp
1705 aval = 0.01558470_dp
1706 ngto(0) = 33
1707 ngto(1) = 24
1708 ngto(2) = 20
1709 ival(1) = 1
1710 CASE (40)
1711 cval = 1.97976190_dp
1712 aval = 0.01705520_dp
1713 ngto(0) = 33
1714 ngto(1) = 24
1715 ngto(2) = 20
1716 ival(1) = 1
1717 CASE (41)
1718 cval = 1.97989290_dp
1719 aval = 0.01527040_dp
1720 ngto(0) = 33
1721 ngto(1) = 24
1722 ngto(2) = 20
1723 ival(1) = 1
1724 CASE (42)
1725 cval = 1.97909240_dp
1726 aval = 0.01879720_dp
1727 ngto(0) = 32
1728 ngto(1) = 24
1729 ngto(2) = 20
1730 ival(1) = 1
1731 CASE (43)
1732 cval = 1.98508430_dp
1733 aval = 0.01497550_dp
1734 ngto(0) = 32
1735 ngto(1) = 24
1736 ngto(2) = 20
1737 ival(1) = 2
1738 ival(2) = 1
1739 CASE (44)
1740 cval = 1.98515010_dp
1741 aval = 0.01856670_dp
1742 ngto(0) = 32
1743 ngto(1) = 24
1744 ngto(2) = 20
1745 ival(1) = 2
1746 ival(2) = 1
1747 CASE (45)
1748 cval = 1.98502970_dp
1749 aval = 0.01487000_dp
1750 ngto(0) = 32
1751 ngto(1) = 24
1752 ngto(2) = 20
1753 ival(1) = 2
1754 ival(2) = 1
1755 CASE (46)
1756 cval = 1.97672850_dp
1757 aval = 0.01762500_dp
1758 ngto(0) = 30
1759 ngto(1) = 24
1760 ngto(2) = 20
1761 ival(0) = 2
1762 ival(1) = 2
1763 ival(2) = 1
1764 CASE (47)
1765 cval = 1.97862730_dp
1766 aval = 0.01863310_dp
1767 ngto(0) = 32
1768 ngto(1) = 24
1769 ngto(2) = 20
1770 ival(1) = 2
1771 ival(2) = 1
1772 CASE (48)
1773 cval = 1.97990020_dp
1774 aval = 0.01347150_dp
1775 ngto(0) = 33
1776 ngto(1) = 24
1777 ngto(2) = 20
1778 ival(1) = 2
1779 ival(2) = 2
1780 CASE (49)
1781 cval = 1.97979410_dp
1782 aval = 0.00890265_dp
1783 ngto(0) = 33
1784 ngto(1) = 27
1785 ngto(2) = 20
1786 ival(0) = 2
1787 ival(2) = 2
1788 CASE (50)
1789 cval = 1.98001000_dp
1790 aval = 0.00895215_dp
1791 ngto(0) = 33
1792 ngto(1) = 27
1793 ngto(2) = 20
1794 ival(0) = 2
1795 ival(2) = 2
1796 CASE (51)
1797 cval = 1.97979980_dp
1798 aval = 0.01490290_dp
1799 ngto(0) = 33
1800 ngto(1) = 26
1801 ngto(2) = 20
1802 ival(1) = 1
1803 ival(2) = 2
1804 CASE (52)
1805 cval = 1.98009310_dp
1806 aval = 0.01490390_dp
1807 ngto(0) = 33
1808 ngto(1) = 26
1809 ngto(2) = 20
1810 ival(1) = 1
1811 ival(2) = 2
1812 CASE (53)
1813 cval = 1.97794750_dp
1814 aval = 0.01425880_dp
1815 ngto(0) = 33
1816 ngto(1) = 26
1817 ngto(2) = 20
1818 ival(0) = 2
1819 ival(1) = 1
1820 ival(2) = 2
1821 CASE (54)
1822 cval = 1.97784450_dp
1823 aval = 0.01430130_dp
1824 ngto(0) = 33
1825 ngto(1) = 26
1826 ngto(2) = 20
1827 ival(0) = 2
1828 ival(1) = 1
1829 ival(2) = 2
1830 CASE (55)
1831 cval = 1.97784450_dp
1832 aval = 0.00499318_dp
1833 ngto(0) = 32
1834 ngto(1) = 25
1835 ngto(2) = 17
1836 ival(0) = 1
1837 ival(1) = 3
1838 ival(2) = 6
1839 CASE (56)
1840 cval = 1.97764820_dp
1841 aval = 0.00500392_dp
1842 ngto(0) = 32
1843 ngto(1) = 25
1844 ngto(2) = 17
1845 ival(0) = 1
1846 ival(1) = 3
1847 ival(2) = 6
1848 CASE (57)
1849 cval = 1.97765150_dp
1850 aval = 0.00557083_dp
1851 ngto(0) = 32
1852 ngto(1) = 25
1853 ngto(2) = 20
1854 ival(0) = 1
1855 ival(1) = 3
1856 ival(2) = 3
1857 CASE (58)
1858 cval = 1.97768750_dp
1859 aval = 0.00547531_dp
1860 ngto(0) = 32
1861 ngto(1) = 25
1862 ngto(2) = 20
1863 ngto(3) = 16
1864 ival(0) = 1
1865 ival(1) = 3
1866 ival(2) = 3
1867 ival(3) = 3
1868 CASE (59)
1869 cval = 1.96986600_dp
1870 aval = 0.00813143_dp
1871 ngto(0) = 32
1872 ngto(1) = 25
1873 ngto(2) = 17
1874 ngto(3) = 16
1875 ival(0) = 1
1876 ival(1) = 3
1877 ival(2) = 6
1878 ival(3) = 4
1879 CASE (60)
1880 cval = 1.97765720_dp
1881 aval = 0.00489201_dp
1882 ngto(0) = 32
1883 ngto(1) = 25
1884 ngto(2) = 17
1885 ngto(3) = 16
1886 ival(0) = 1
1887 ival(1) = 3
1888 ival(2) = 6
1889 ival(3) = 4
1890 CASE (61)
1891 cval = 1.97768120_dp
1892 aval = 0.00499000_dp
1893 ngto(0) = 32
1894 ngto(1) = 25
1895 ngto(2) = 17
1896 ngto(3) = 16
1897 ival(0) = 1
1898 ival(1) = 3
1899 ival(2) = 6
1900 ival(3) = 4
1901 CASE (62)
1902 cval = 1.97745700_dp
1903 aval = 0.00615587_dp
1904 ngto(0) = 32
1905 ngto(1) = 25
1906 ngto(2) = 17
1907 ngto(3) = 16
1908 ival(0) = 1
1909 ival(1) = 3
1910 ival(2) = 6
1911 ival(3) = 4
1912 CASE (63)
1913 cval = 1.97570240_dp
1914 aval = 0.00769959_dp
1915 ngto(0) = 32
1916 ngto(1) = 25
1917 ngto(2) = 17
1918 ngto(3) = 16
1919 ival(0) = 1
1920 ival(1) = 3
1921 ival(2) = 6
1922 ival(3) = 4
1923 CASE (64)
1924 cval = 1.97629350_dp
1925 aval = 0.00706610_dp
1926 ngto(0) = 32
1927 ngto(1) = 25
1928 ngto(2) = 20
1929 ngto(3) = 16
1930 ival(0) = 1
1931 ival(1) = 3
1932 ival(2) = 3
1933 ival(3) = 4
1934 CASE (65)
1935 cval = 1.96900000_dp
1936 aval = 0.01019150_dp
1937 ngto(0) = 32
1938 ngto(1) = 26
1939 ngto(2) = 18
1940 ngto(3) = 16
1941 ival(0) = 1
1942 ival(1) = 3
1943 ival(2) = 6
1944 ival(3) = 4
1945 CASE (66)
1946 cval = 1.97350000_dp
1947 aval = 0.01334320_dp
1948 ngto(0) = 33
1949 ngto(1) = 26
1950 ngto(2) = 18
1951 ngto(3) = 16
1952 ival(0) = 1
1953 ival(1) = 3
1954 ival(2) = 6
1955 ival(3) = 4
1956 CASE (67)
1957 cval = 1.97493000_dp
1958 aval = 0.01331360_dp
1959 ngto(0) = 32
1960 ngto(1) = 24
1961 ngto(2) = 17
1962 ngto(3) = 14
1963 ival(1) = 2
1964 ival(2) = 5
1965 ival(3) = 4
1966 CASE (68)
1967 cval = 1.97597670_dp
1968 aval = 0.01434040_dp
1969 ngto(0) = 32
1970 ngto(1) = 24
1971 ngto(2) = 17
1972 ngto(3) = 14
1973 ival(0) = 0
1974 ival(1) = 2
1975 ival(2) = 5
1976 ival(3) = 4
1977 CASE (69)
1978 cval = 1.97809240_dp
1979 aval = 0.01529430_dp
1980 ngto(0) = 32
1981 ngto(1) = 24
1982 ngto(2) = 17
1983 ngto(3) = 14
1984 ival(0) = 0
1985 ival(1) = 2
1986 ival(2) = 5
1987 ival(3) = 4
1988 CASE (70)
1989 cval = 1.97644360_dp
1990 aval = 0.01312770_dp
1991 ngto(0) = 32
1992 ngto(1) = 24
1993 ngto(2) = 17
1994 ngto(3) = 14
1995 ival(0) = 0
1996 ival(1) = 2
1997 ival(2) = 5
1998 ival(3) = 4
1999 CASE (71)
2000 cval = 1.96998000_dp
2001 aval = 0.01745150_dp
2002 ngto(0) = 31
2003 ngto(1) = 24
2004 ngto(2) = 20
2005 ngto(3) = 14
2006 ival(0) = 1
2007 ival(1) = 2
2008 ival(2) = 2
2009 ival(3) = 4
2010 CASE (72)
2011 cval = 1.97223830_dp
2012 aval = 0.01639750_dp
2013 ngto(0) = 31
2014 ngto(1) = 24
2015 ngto(2) = 20
2016 ngto(3) = 14
2017 ival(0) = 1
2018 ival(1) = 2
2019 ival(2) = 2
2020 ival(3) = 4
2021 CASE (73)
2022 cval = 1.97462110_dp
2023 aval = 0.01603680_dp
2024 ngto(0) = 31
2025 ngto(1) = 24
2026 ngto(2) = 20
2027 ngto(3) = 14
2028 ival(0) = 1
2029 ival(1) = 2
2030 ival(2) = 2
2031 ival(3) = 4
2032 CASE (74)
2033 cval = 1.97756000_dp
2034 aval = 0.02030570_dp
2035 ngto(0) = 31
2036 ngto(1) = 24
2037 ngto(2) = 20
2038 ngto(3) = 14
2039 ival(0) = 1
2040 ival(1) = 2
2041 ival(2) = 2
2042 ival(3) = 4
2043 CASE (75)
2044 cval = 1.97645760_dp
2045 aval = 0.02057180_dp
2046 ngto(0) = 31
2047 ngto(1) = 24
2048 ngto(2) = 20
2049 ngto(3) = 14
2050 ival(0) = 1
2051 ival(1) = 2
2052 ival(2) = 2
2053 ival(3) = 4
2054 CASE (76)
2055 cval = 1.97725820_dp
2056 aval = 0.02058210_dp
2057 ngto(0) = 32
2058 ngto(1) = 24
2059 ngto(2) = 20
2060 ngto(3) = 15
2061 ival(0) = 0
2062 ival(1) = 2
2063 ival(2) = 2
2064 ival(3) = 4
2065 CASE (77)
2066 cval = 1.97749380_dp
2067 aval = 0.02219380_dp
2068 ngto(0) = 32
2069 ngto(1) = 24
2070 ngto(2) = 20
2071 ngto(3) = 15
2072 ival(0) = 0
2073 ival(1) = 2
2074 ival(2) = 2
2075 ival(3) = 4
2076 CASE (78)
2077 cval = 1.97946280_dp
2078 aval = 0.02216280_dp
2079 ngto(0) = 32
2080 ngto(1) = 24
2081 ngto(2) = 20
2082 ngto(3) = 15
2083 ival(0) = 0
2084 ival(1) = 2
2085 ival(2) = 2
2086 ival(3) = 4
2087 CASE (79)
2088 cval = 1.97852130_dp
2089 aval = 0.02168500_dp
2090 ngto(0) = 32
2091 ngto(1) = 24
2092 ngto(2) = 20
2093 ngto(3) = 15
2094 ival(0) = 0
2095 ival(1) = 2
2096 ival(2) = 2
2097 ival(3) = 4
2098 CASE (80)
2099 cval = 1.98045190_dp
2100 aval = 0.02177860_dp
2101 ngto(0) = 32
2102 ngto(1) = 24
2103 ngto(2) = 20
2104 ngto(3) = 15
2105 ival(0) = 0
2106 ival(1) = 2
2107 ival(2) = 2
2108 ival(3) = 4
2109 CASE (81)
2110 cval = 1.97000000_dp
2111 aval = 0.02275000_dp
2112 ngto(0) = 31
2113 ngto(1) = 25
2114 ngto(2) = 18
2115 ngto(3) = 13
2116 ival(0) = 1
2117 ival(1) = 0
2118 ival(2) = 3
2119 ival(3) = 6
2120 CASE (82)
2121 cval = 1.97713580_dp
2122 aval = 0.02317030_dp
2123 ngto(0) = 31
2124 ngto(1) = 27
2125 ngto(2) = 18
2126 ngto(3) = 13
2127 ival(0) = 1
2128 ival(1) = 0
2129 ival(2) = 3
2130 ival(3) = 6
2131 CASE (83)
2132 cval = 1.97537880_dp
2133 aval = 0.02672860_dp
2134 ngto(0) = 32
2135 ngto(1) = 27
2136 ngto(2) = 17
2137 ngto(3) = 13
2138 ival(0) = 1
2139 ival(1) = 0
2140 ival(2) = 3
2141 ival(3) = 6
2142 CASE (84)
2143 cval = 1.97545360_dp
2144 aval = 0.02745360_dp
2145 ngto(0) = 31
2146 ngto(1) = 27
2147 ngto(2) = 17
2148 ngto(3) = 13
2149 ival(0) = 1
2150 ival(1) = 0
2151 ival(2) = 3
2152 ival(3) = 6
2153 CASE (85)
2154 cval = 1.97338370_dp
2155 aval = 0.02616310_dp
2156 ngto(0) = 32
2157 ngto(1) = 27
2158 ngto(2) = 19
2159 ngto(3) = 13
2160 ival(0) = 1
2161 ival(1) = 0
2162 ival(2) = 3
2163 ival(3) = 6
2164 CASE (86)
2165 cval = 1.97294240_dp
2166 aval = 0.02429220_dp
2167 ngto(0) = 32
2168 ngto(1) = 27
2169 ngto(2) = 19
2170 ngto(3) = 13
2171 ival(0) = 1
2172 ival(1) = 0
2173 ival(2) = 3
2174 ival(3) = 6
2175 CASE (87:106) ! these numbers are an educated guess
2176 cval = 1.98000000_dp
2177 aval = 0.01400000_dp
2178 ngto(0) = 34
2179 ngto(1) = 28
2180 ngto(2) = 20
2181 ngto(3) = 15
2182 ival(0) = 0
2183 ival(1) = 0
2184 ival(2) = 3
2185 ival(3) = 6
2186 CASE DEFAULT
2187 cpabort("No geometrical basis set data are available for the selected atom number.")
2188 END SELECT
2189
2190 END SUBROUTINE clementi_geobas
2191
2192! **************************************************************************************************
2193!> \brief ...
2194!> \param element_symbol ...
2195!> \param basis ...
2196!> \param basis_set_name ...
2197!> \param basis_set_file ...
2198!> \param basis_section ...
2199! **************************************************************************************************
2200 SUBROUTINE read_basis_set(element_symbol, basis, basis_set_name, basis_set_file, &
2201 basis_section)
2202
2203 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
2204 TYPE(atom_basis_type), INTENT(INOUT) :: basis
2205 CHARACTER(LEN=*), INTENT(IN) :: basis_set_name, basis_set_file
2206 TYPE(section_vals_type), POINTER :: basis_section
2207
2208 INTEGER, PARAMETER :: maxpri = 40, maxset = 20
2209
2210 CHARACTER(len=20*default_string_length) :: line_att
2211 CHARACTER(LEN=240) :: line
2212 CHARACTER(LEN=242) :: line2
2213 CHARACTER(LEN=LEN(basis_set_name)) :: bsname
2214 CHARACTER(LEN=LEN(basis_set_name)+2) :: bsname2
2215 CHARACTER(LEN=LEN(element_symbol)) :: symbol
2216 CHARACTER(LEN=LEN(element_symbol)+2) :: symbol2
2217 INTEGER :: i, ii, ipgf, irep, iset, ishell, j, k, &
2218 lshell, nj, nmin, ns, nset, strlen1, &
2219 strlen2
2220 INTEGER, DIMENSION(maxpri, maxset) :: l
2221 INTEGER, DIMENSION(maxset) :: lmax, lmin, n, npgf, nshell
2222 LOGICAL :: found, is_ok, match, read_from_input
2223 REAL(dp) :: expzet, gcca, prefac, zeta
2224 REAL(dp), DIMENSION(maxpri, maxpri, maxset) :: gcc
2225 REAL(dp), DIMENSION(maxpri, maxset) :: zet
2226 TYPE(cp_sll_val_type), POINTER :: list
2227 TYPE(val_type), POINTER :: val
2228
2229 bsname = basis_set_name
2230 symbol = element_symbol
2231 irep = 0
2232
2233 nset = 0
2234 lmin = 0
2235 lmax = 0
2236 npgf = 0
2237 n = 0
2238 l = 0
2239 zet = 0._dp
2240 gcc = 0._dp
2241
2242 read_from_input = .false.
2243 CALL section_vals_get(basis_section, explicit=read_from_input)
2244 IF (read_from_input) THEN
2245 NULLIFY (list, val)
2246 CALL section_vals_list_get(basis_section, "_DEFAULT_KEYWORD_", list=list)
2247 CALL uppercase(symbol)
2248 CALL uppercase(bsname)
2249 is_ok = cp_sll_val_next(list, val)
2250 cpassert(is_ok)
2251 CALL val_get(val, c_val=line_att)
2252 READ (line_att, *) nset
2253 cpassert(nset <= maxset)
2254 DO iset = 1, nset
2255 is_ok = cp_sll_val_next(list, val)
2256 cpassert(is_ok)
2257 CALL val_get(val, c_val=line_att)
2258 READ (line_att, *) n(iset)
2259 CALL remove_word(line_att)
2260 READ (line_att, *) lmin(iset)
2261 CALL remove_word(line_att)
2262 READ (line_att, *) lmax(iset)
2263 CALL remove_word(line_att)
2264 READ (line_att, *) npgf(iset)
2265 CALL remove_word(line_att)
2266 cpassert(npgf(iset) <= maxpri)
2267 nshell(iset) = 0
2268 DO lshell = lmin(iset), lmax(iset)
2269 nmin = n(iset) + lshell - lmin(iset)
2270 READ (line_att, *) ishell
2271 CALL remove_word(line_att)
2272 nshell(iset) = nshell(iset) + ishell
2273 DO i = 1, ishell
2274 l(nshell(iset) - ishell + i, iset) = lshell
2275 END DO
2276 END DO
2277 cpassert(len_trim(line_att) == 0)
2278 DO ipgf = 1, npgf(iset)
2279 is_ok = cp_sll_val_next(list, val)
2280 cpassert(is_ok)
2281 CALL val_get(val, c_val=line_att)
2282 READ (line_att, *) zet(ipgf, iset), (gcc(ipgf, ishell, iset), ishell=1, nshell(iset))
2283 END DO
2284 END DO
2285 ELSE
2286 block
2287 TYPE(cp_parser_type) :: parser
2288 CALL parser_create(parser, basis_set_file)
2289 ! Search for the requested basis set in the basis set file
2290 ! until the basis set is found or the end of file is reached
2291 search_loop: DO
2292 CALL parser_search_string(parser, trim(bsname), .true., found, line)
2293 IF (found) THEN
2294 CALL uppercase(symbol)
2295 CALL uppercase(bsname)
2296 match = .false.
2297 CALL uppercase(line)
2298 ! Check both the element symbol and the basis set name
2299 line2 = " "//line//" "
2300 symbol2 = " "//trim(symbol)//" "
2301 bsname2 = " "//trim(bsname)//" "
2302 strlen1 = len_trim(symbol2) + 1
2303 strlen2 = len_trim(bsname2) + 1
2304
2305 IF ((index(line2, symbol2(:strlen1)) > 0) .AND. &
2306 (index(line2, bsname2(:strlen2)) > 0)) match = .true.
2307
2308 IF (match) THEN
2309 ! Read the basis set information
2310 CALL parser_get_object(parser, nset, newline=.true.)
2311 cpassert(nset <= maxset)
2312 DO iset = 1, nset
2313 CALL parser_get_object(parser, n(iset), newline=.true.)
2314 CALL parser_get_object(parser, lmin(iset))
2315 CALL parser_get_object(parser, lmax(iset))
2316 CALL parser_get_object(parser, npgf(iset))
2317 cpassert(npgf(iset) <= maxpri)
2318 nshell(iset) = 0
2319 DO lshell = lmin(iset), lmax(iset)
2320 nmin = n(iset) + lshell - lmin(iset)
2321 CALL parser_get_object(parser, ishell)
2322 nshell(iset) = nshell(iset) + ishell
2323 DO i = 1, ishell
2324 l(nshell(iset) - ishell + i, iset) = lshell
2325 END DO
2326 END DO
2327 DO ipgf = 1, npgf(iset)
2328 CALL parser_get_object(parser, zet(ipgf, iset), newline=.true.)
2329 DO ishell = 1, nshell(iset)
2330 CALL parser_get_object(parser, gcc(ipgf, ishell, iset))
2331 END DO
2332 END DO
2333 END DO
2334
2335 EXIT search_loop
2336
2337 END IF
2338 ELSE
2339 ! Stop program, if the end of file is reached
2340 cpabort("End of file reached and the requested basis set was not found.")
2341 END IF
2342
2343 END DO search_loop
2344
2345 CALL parser_release(parser)
2346 END block
2347 END IF
2348
2349 ! fill in the basis data structures
2350 basis%nprim = 0
2351 basis%nbas = 0
2352 DO i = 1, nset
2353 DO j = lmin(i), min(lmax(i), lmat)
2354 basis%nprim(j) = basis%nprim(j) + npgf(i)
2355 END DO
2356 DO j = 1, nshell(i)
2357 k = l(j, i)
2358 IF (k <= lmat) basis%nbas(k) = basis%nbas(k) + 1
2359 END DO
2360 END DO
2361
2362 nj = maxval(basis%nprim)
2363 ns = maxval(basis%nbas)
2364 ALLOCATE (basis%am(nj, 0:lmat))
2365 basis%am = 0._dp
2366 ALLOCATE (basis%cm(nj, ns, 0:lmat))
2367 basis%cm = 0._dp
2368
2369 DO j = 0, lmat
2370 nj = 0
2371 ns = 0
2372 DO i = 1, nset
2373 IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
2374 DO ipgf = 1, npgf(i)
2375 basis%am(nj + ipgf, j) = zet(ipgf, i)
2376 END DO
2377 DO ii = 1, nshell(i)
2378 IF (l(ii, i) == j) THEN
2379 ns = ns + 1
2380 DO ipgf = 1, npgf(i)
2381 basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
2382 END DO
2383 END IF
2384 END DO
2385 nj = nj + npgf(i)
2386 END IF
2387 END DO
2388 END DO
2389
2390 ! Normalization
2391 DO j = 0, lmat
2392 expzet = 0.25_dp*real(2*j + 3, dp)
2393 prefac = sqrt(rootpi/2._dp**(j + 2)*dfac(2*j + 1))
2394 DO ipgf = 1, basis%nprim(j)
2395 DO ii = 1, basis%nbas(j)
2396 gcca = basis%cm(ipgf, ii, j)
2397 zeta = 2._dp*basis%am(ipgf, j)
2398 basis%cm(ipgf, ii, j) = zeta**expzet*gcca/prefac
2399 END DO
2400 END DO
2401 END DO
2402
2403 END SUBROUTINE read_basis_set
2404
2405! **************************************************************************************************
2406!> \brief ...
2407!> \param optimization ...
2408!> \param opt_section ...
2409! **************************************************************************************************
2410 SUBROUTINE read_atom_opt_section(optimization, opt_section)
2411 TYPE(atom_optimization_type), INTENT(INOUT) :: optimization
2412 TYPE(section_vals_type), POINTER :: opt_section
2413
2414 INTEGER :: miter, ndiis
2415 REAL(kind=dp) :: damp, eps_diis, eps_scf
2416
2417 CALL section_vals_val_get(opt_section, "MAX_ITER", i_val=miter)
2418 CALL section_vals_val_get(opt_section, "EPS_SCF", r_val=eps_scf)
2419 CALL section_vals_val_get(opt_section, "N_DIIS", i_val=ndiis)
2420 CALL section_vals_val_get(opt_section, "EPS_DIIS", r_val=eps_diis)
2421 CALL section_vals_val_get(opt_section, "DAMPING", r_val=damp)
2422
2423 optimization%max_iter = miter
2424 optimization%eps_scf = eps_scf
2425 optimization%n_diis = ndiis
2426 optimization%eps_diis = eps_diis
2427 optimization%damping = damp
2428
2429 END SUBROUTINE read_atom_opt_section
2430! **************************************************************************************************
2431!> \brief ...
2432!> \param potential ...
2433!> \param potential_section ...
2434!> \param zval ...
2435! **************************************************************************************************
2436 SUBROUTINE init_atom_potential(potential, potential_section, zval)
2437 TYPE(atom_potential_type), INTENT(INOUT) :: potential
2438 TYPE(section_vals_type), POINTER :: potential_section
2439 INTEGER, INTENT(IN) :: zval
2440
2441 CHARACTER(LEN=default_string_length) :: pseudo_fn, pseudo_name
2442 INTEGER :: ic
2443 REAL(dp), DIMENSION(:), POINTER :: convals
2444 TYPE(section_vals_type), POINTER :: ecp_potential_section, &
2445 gth_potential_section
2446
2447 IF (zval > 0) THEN
2448 CALL section_vals_val_get(potential_section, "PSEUDO_TYPE", i_val=potential%ppot_type)
2449
2450 SELECT CASE (potential%ppot_type)
2451 CASE (gth_pseudo)
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)
2457 CASE (ecp_pseudo)
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")
2461 CALL read_ecp_potential(ptable(zval)%symbol, potential%ecp_pot, &
2462 pseudo_name, pseudo_fn, ecp_potential_section)
2463 CASE (upf_pseudo)
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
2468 CASE (sgp_pseudo)
2469 cpabort("Pseudopotential type SGP is not implemented.")
2470 CASE (no_pseudo)
2471 ! do nothing
2472 CASE DEFAULT
2473 cpabort("Invalid pseudopotential type selected. Check the code!")
2474 END SELECT
2475 ELSE
2476 potential%ppot_type = no_pseudo
2477 END IF
2478
2479 ! confinement
2480 NULLIFY (convals)
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)
2496 ELSE
2497 potential%rcon = 4.0_dp
2498 END IF
2499 IF (SIZE(convals) >= 3) THEN
2500 potential%scon = convals(3)
2501 ELSE
2502 potential%scon = 2.0_dp
2503 END IF
2504 ELSE
2505 potential%confinement = .false.
2506 END IF
2507 ELSE
2508 potential%confinement = .false.
2509 END IF
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)
2521 END IF
2522 IF (SIZE(convals) >= 3) THEN
2523 potential%scon = convals(3)
2524 END IF
2525 ELSE
2526 potential%confinement = .false.
2527 END IF
2528 END IF
2529 END IF
2530
2531 END SUBROUTINE init_atom_potential
2532! **************************************************************************************************
2533!> \brief ...
2534!> \param potential ...
2535! **************************************************************************************************
2536 SUBROUTINE release_atom_potential(potential)
2537 TYPE(atom_potential_type), INTENT(INOUT) :: potential
2538
2539 potential%confinement = .false.
2540
2541 CALL atom_release_upf(potential%upf_pot)
2542
2543 END SUBROUTINE release_atom_potential
2544! **************************************************************************************************
2545!> \brief ...
2546!> \param element_symbol ...
2547!> \param potential ...
2548!> \param pseudo_name ...
2549!> \param pseudo_file ...
2550!> \param potential_section ...
2551! **************************************************************************************************
2552 SUBROUTINE read_gth_potential(element_symbol, potential, pseudo_name, pseudo_file, &
2553 potential_section)
2554
2555 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
2556 TYPE(atom_gthpot_type), INTENT(INOUT) :: potential
2557 CHARACTER(LEN=*), INTENT(IN) :: pseudo_name, pseudo_file
2558 TYPE(section_vals_type), POINTER :: potential_section
2559
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, &
2568 strlen2
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
2573
2574 elec_conf = 0
2575
2576 apname = pseudo_name
2577 symbol = element_symbol
2578
2579 potential%symbol = symbol
2580 potential%pname = apname
2581 potential%econf = 0
2582 potential%rc = 0._dp
2583 potential%ncl = 0
2584 potential%cl = 0._dp
2585 potential%nl = 0
2586 potential%rcnl = 0._dp
2587 potential%hnl = 0._dp
2588 potential%soc = .false.
2589 potential%knl = 0._dp
2590
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
2597
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)
2604 ! Read the electronic configuration, not used here
2605 l = 0
2606 is_ok = cp_sll_val_next(list, val)
2607 cpassert(is_ok)
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)
2612 l = l + 1
2613 READ (line_att, *) elec_conf(l)
2614 CALL remove_word(line_att)
2615 END DO
2616 potential%econf(0:lmat) = elec_conf(0:lmat)
2617 potential%zion = real(sum(elec_conf), dp)
2618 ! Read r(loc) to define the exponent of the core charge
2619 is_ok = cp_sll_val_next(list, val)
2620 cpassert(is_ok)
2621 CALL val_get(val, c_val=line_att)
2622 READ (line_att, *) potential%rc
2623 CALL remove_word(line_att)
2624 ! Read the parameters for the local part of the GTH pseudopotential (ppl)
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)
2630 END DO
2631 ! Check for the next entry: LPOT, NLCC, LSD, or ppnl
2632 DO
2633 is_ok = cp_sll_val_next(list, val)
2634 cpassert(is_ok)
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)
2642 cpassert(is_ok)
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)
2651 END DO
2652 END DO
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)
2659 cpassert(is_ok)
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)
2667 !make cp2k compatible with bigdft
2668 potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2669 CALL remove_word(line_att)
2670 END DO
2671 END DO
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)
2678 cpassert(is_ok)
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)
2687 END DO
2688 END DO
2689 ELSE
2690 EXIT
2691 END IF
2692 END DO
2693 ! Read the parameters for the non-local part of the GTH pseudopotential (ppnl)
2694 READ (line_att, *) nlmax
2695 CALL remove_word(line_att)
2696 IF (index(line_att, "SOC") /= 0) potential%soc = .true.
2697 IF (nlmax > 0) THEN
2698 ! Load the parameter for nlmax non-local projectors
2699 DO l = 0, nlmax - 1
2700 is_ok = cp_sll_val_next(list, val)
2701 cpassert(is_ok)
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)
2708 IF (i == 1) THEN
2709 READ (line_att, *) potential%hnl(1, 1, l)
2710 CALL remove_word(line_att)
2711 ELSE
2712 cpassert(len_trim(line_att) == 0)
2713 is_ok = cp_sll_val_next(list, val)
2714 cpassert(is_ok)
2715 CALL val_get(val, c_val=line_att)
2716 READ (line_att, *) potential%hnl(i, i, l)
2717 CALL remove_word(line_att)
2718 END IF
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)
2723 END DO
2724 END DO
2725 IF (potential%soc .AND. l /= 0) THEN
2726 is_ok = cp_sll_val_next(list, val)
2727 cpassert(is_ok)
2728 CALL val_get(val, c_val=line_att)
2729 DO i = 1, potential%nl(l)
2730 IF (i == 1) THEN
2731 READ (line_att, *) potential%knl(1, 1, l)
2732 CALL remove_word(line_att)
2733 ELSE
2734 cpassert(len_trim(line_att) == 0)
2735 is_ok = cp_sll_val_next(list, val)
2736 cpassert(is_ok)
2737 CALL val_get(val, c_val=line_att)
2738 READ (line_att, *) potential%knl(i, i, l)
2739 CALL remove_word(line_att)
2740 END IF
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)
2745 END DO
2746 END DO
2747 END IF
2748 cpassert(len_trim(line_att) == 0)
2749 END DO
2750 END IF
2751 ELSE
2752 block
2753 TYPE(cp_parser_type) :: parser
2754 CALL parser_create(parser, pseudo_file)
2755
2756 search_loop: DO
2757 CALL parser_search_string(parser, trim(apname), .true., found, line)
2758 IF (found) THEN
2759 CALL uppercase(symbol)
2760 CALL uppercase(apname)
2761 ! Check both the element symbol and the atomic potential name
2762 match = .false.
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
2769
2770 IF ((index(line2, symbol2(:strlen1)) > 0) .AND. &
2771 (index(line2, apname2(:strlen2)) > 0)) match = .true.
2772
2773 IF (match) THEN
2774 ! Read the electronic configuration
2775 l = 0
2776 CALL parser_get_object(parser, elec_conf(l), newline=.true.)
2777 DO WHILE (parser_test_next_token(parser) == "INT")
2778 l = l + 1
2779 CALL parser_get_object(parser, elec_conf(l))
2780 END DO
2781 potential%econf(0:lmat) = elec_conf(0:lmat)
2782 potential%zion = real(sum(elec_conf), dp)
2783 ! Read r(loc) to define the exponent of the core charge
2784 CALL parser_get_object(parser, potential%rc, newline=.true.)
2785 ! Read the parameters for the local part of the GTH pseudopotential (ppl)
2786 CALL parser_get_object(parser, potential%ncl)
2787 DO i = 1, potential%ncl
2788 CALL parser_get_object(parser, potential%cl(i))
2789 END DO
2790 ! Extended type input
2791 DO
2792 CALL parser_get_next_line(parser, 1)
2793 IF (parser_test_next_token(parser) == "INT") THEN
2794 EXIT
2795 ELSE IF (parser_test_next_token(parser) == "STR") THEN
2796 CALL parser_get_object(parser, line)
2797 IF (index(line, "LPOT") /= 0) THEN
2798 ! local potential
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))
2806 END DO
2807 END DO
2808 ELSE IF (index(line, "NLCC") /= 0) THEN
2809 ! NLCC
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))
2817 !make cp2k compatible with bigdft
2818 potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2819 END DO
2820 END DO
2821 ELSE IF (index(line, "LSD") /= 0) THEN
2822 ! LSD potential
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))
2830 END DO
2831 END DO
2832 ELSE
2833 cpabort("Parsing of extended potential type failed.")
2834 END IF
2835 ELSE
2836 cpabort("Invalid input token found.")
2837 END IF
2838 END DO
2839 ! Read the parameters for the non-local part of the GTH pseudopotential (ppnl)
2840 CALL parser_get_object(parser, nlmax)
2841 IF (nlmax > 0) THEN
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.
2845 END IF
2846 ! Load the parameter for n non-local projectors
2847 DO l = 0, nlmax - 1
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)
2851 IF (i == 1) THEN
2852 CALL parser_get_object(parser, potential%hnl(i, i, l))
2853 ELSE
2854 CALL parser_get_object(parser, potential%hnl(i, i, l), newline=.true.)
2855 END IF
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)
2859 END DO
2860 END DO
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)
2867 END DO
2868 END DO
2869 END IF
2870 END DO
2871 END IF
2872 EXIT search_loop
2873 END IF
2874 ELSE
2875 ! Stop program, if the end of file is reached
2876 cpabort("End of file reached unexpectedly")
2877 END IF
2878
2879 END DO search_loop
2880
2881 CALL parser_release(parser)
2882 END block
2883 END IF
2884
2885 END SUBROUTINE read_gth_potential
2886! **************************************************************************************************
2887!> \brief ...
2888!> \param element_symbol ...
2889!> \param potential ...
2890!> \param pseudo_name ...
2891!> \param pseudo_file ...
2892!> \param potential_section ...
2893! **************************************************************************************************
2894 SUBROUTINE read_ecp_potential(element_symbol, potential, pseudo_name, pseudo_file, &
2895 potential_section)
2896
2897 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
2898 TYPE(atom_ecppot_type), INTENT(INOUT) :: potential
2899 CHARACTER(LEN=*), INTENT(IN) :: pseudo_name, pseudo_file
2900 TYPE(section_vals_type), POINTER :: potential_section
2901
2902 CHARACTER(LEN=240) :: line
2903 CHARACTER(len=5*default_string_length) :: line_att
2904 CHARACTER(LEN=LEN(element_symbol)+1) :: symbol
2905 CHARACTER(LEN=LEN(pseudo_name)) :: apname
2906 INTEGER :: i, ic, l, ncore, nel
2907 LOGICAL :: found, is_ok, read_from_input
2908 TYPE(cp_sll_val_type), POINTER :: list
2909 TYPE(val_type), POINTER :: val
2910
2911 apname = pseudo_name
2912 symbol = element_symbol
2913 CALL get_ptable_info(symbol, number=ncore)
2914
2915 potential%symbol = symbol
2916 potential%pname = apname
2917 potential%econf = 0
2918 potential%zion = 0
2919 potential%lmax = -1
2920 potential%nloc = 0
2921 potential%nrloc = 0
2922 potential%aloc = 0.0_dp
2923 potential%bloc = 0.0_dp
2924 potential%npot = 0
2925 potential%nrpot = 0
2926 potential%apot = 0.0_dp
2927 potential%bpot = 0.0_dp
2928
2929 read_from_input = .false.
2930 CALL section_vals_get(potential_section, explicit=read_from_input)
2931 IF (read_from_input) THEN
2932 CALL section_vals_list_get(potential_section, "_DEFAULT_KEYWORD_", list=list)
2933 ! number of electrons (mandatory line)
2934 is_ok = cp_sll_val_next(list, val)
2935 cpassert(is_ok)
2936 CALL val_get(val, c_val=line_att)
2937 CALL remove_word(line_att)
2938 CALL remove_word(line_att)
2939 ! read number of electrons
2940 READ (line_att, *) nel
2941 potential%zion = real(ncore - nel, kind=dp)
2942 ! local potential (mandatory block)
2943 is_ok = cp_sll_val_next(list, val)
2944 cpassert(is_ok)
2945 CALL val_get(val, c_val=line_att)
2946 DO i = 1, 10
2947 IF (.NOT. cp_sll_val_next(list, val)) EXIT
2948 CALL val_get(val, c_val=line_att)
2949 IF (index(line_att, element_symbol) == 0) THEN
2950 potential%nloc = potential%nloc + 1
2951 ic = potential%nloc
2952 READ (line_att, *) potential%nrloc(ic), potential%bloc(ic), potential%aloc(ic)
2953 ELSE
2954 EXIT
2955 END IF
2956 END DO
2957 ! read potentials
2958 DO
2959 CALL val_get(val, c_val=line_att)
2960 IF (index(line_att, element_symbol) == 0) THEN
2961 potential%npot(l) = potential%npot(l) + 1
2962 ic = potential%npot(l)
2963 READ (line_att, *) potential%nrpot(ic, l), potential%bpot(ic, l), potential%apot(ic, l)
2964 ELSE
2965 potential%lmax = potential%lmax + 1
2966 l = potential%lmax
2967 END IF
2968 IF (.NOT. cp_sll_val_next(list, val)) EXIT
2969 END DO
2970
2971 ELSE
2972 block
2973 TYPE(cp_parser_type) :: parser
2974 CALL parser_create(parser, pseudo_file)
2975
2976 search_loop: DO
2977 CALL parser_search_string(parser, trim(apname), .true., found, line)
2978 IF (found) THEN
2979 match_loop: DO
2980 CALL parser_get_object(parser, line, newline=.true.)
2981 IF (trim(line) == element_symbol) THEN
2982 CALL parser_get_object(parser, line, lower_to_upper=.true.)
2983 cpassert(trim(line) == "NELEC")
2984 ! read number of electrons
2985 CALL parser_get_object(parser, nel)
2986 potential%zion = real(ncore - nel, kind=dp)
2987 ! read local potential flag line "<XX> ul"
2988 CALL parser_get_object(parser, line, newline=.true.)
2989 ! read local potential
2990 DO i = 1, 15
2991 CALL parser_read_line(parser, 1)
2992 IF (parser_test_next_token(parser) == "STR") EXIT
2993 potential%nloc = potential%nloc + 1
2994 ic = potential%nloc
2995 CALL parser_get_object(parser, potential%nrloc(ic))
2996 CALL parser_get_object(parser, potential%bloc(ic))
2997 CALL parser_get_object(parser, potential%aloc(ic))
2998 END DO
2999 ! read potentials (start with l loop)
3000 DO l = 0, 15
3001 CALL parser_get_object(parser, symbol)
3002 IF (symbol == element_symbol) THEN
3003 ! new l block
3004 potential%lmax = potential%lmax + 1
3005 DO i = 1, 15
3006 CALL parser_read_line(parser, 1)
3007 IF (parser_test_next_token(parser) == "STR") EXIT
3008 potential%npot(l) = potential%npot(l) + 1
3009 ic = potential%npot(l)
3010 CALL parser_get_object(parser, potential%nrpot(ic, l))
3011 CALL parser_get_object(parser, potential%bpot(ic, l))
3012 CALL parser_get_object(parser, potential%apot(ic, l))
3013 END DO
3014 ELSE
3015 EXIT
3016 END IF
3017 END DO
3018 EXIT search_loop
3019 ELSE IF (line == "END") THEN
3020 cpabort("Element not found in ECP library")
3021 END IF
3022 END DO match_loop
3023 ELSE
3024 cpabort("ECP type not found in library")
3025 END IF
3026
3027 END DO search_loop
3028
3029 CALL parser_release(parser)
3030 END block
3031 END IF
3032
3033 ! set up econf
3034 potential%econf(0:3) = ptable(ncore)%e_conv(0:3)
3035 SELECT CASE (nel)
3036 CASE DEFAULT
3037 cpabort("Unknown Core State")
3038 CASE (0)
3039 CASE (2)
3040 potential%econf(0:3) = potential%econf(0:3) - ptable(2)%e_conv(0:3)
3041 CASE (10)
3042 potential%econf(0:3) = potential%econf(0:3) - ptable(10)%e_conv(0:3)
3043 CASE (18)
3044 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
3045 CASE (28)
3046 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
3047 potential%econf(2) = potential%econf(2) - 10
3048 CASE (36)
3049 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3050 CASE (46)
3051 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3052 potential%econf(2) = potential%econf(2) - 10
3053 CASE (54)
3054 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3055 CASE (60)
3056 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3057 potential%econf(2) = potential%econf(2) - 10
3058 potential%econf(3) = potential%econf(3) - 14
3059 CASE (68)
3060 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3061 potential%econf(3) = potential%econf(3) - 14
3062 CASE (78)
3063 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3064 potential%econf(2) = potential%econf(2) - 10
3065 potential%econf(3) = potential%econf(3) - 14
3066 END SELECT
3067 !
3068 cpassert(all(potential%econf >= 0))
3069
3070 END SUBROUTINE read_ecp_potential
3071! **************************************************************************************************
3072!> \brief ...
3073!> \param grid1 ...
3074!> \param grid2 ...
3075!> \return ...
3076! **************************************************************************************************
3077 FUNCTION atom_compare_grids(grid1, grid2) RESULT(is_equal)
3078 TYPE(grid_atom_type) :: grid1, grid2
3079 LOGICAL :: is_equal
3080
3081 INTEGER :: i
3082 REAL(kind=dp) :: dr, dw
3083
3084 is_equal = .true.
3085 IF (grid1%nr == grid2%nr) THEN
3086 DO i = 1, grid2%nr
3087 dr = abs(grid1%rad(i) - grid2%rad(i))
3088 dw = abs(grid1%wr(i) - grid2%wr(i))
3089 IF (dr + dw > 1.0e-12_dp) THEN
3090 is_equal = .false.
3091 EXIT
3092 END IF
3093 END DO
3094 ELSE
3095 is_equal = .false.
3096 END IF
3097
3098 END FUNCTION atom_compare_grids
3099! **************************************************************************************************
3100
3101END MODULE atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
integer, parameter, public num_basis
Definition atom_types.F:69
subroutine, public read_atom_opt_section(optimization, opt_section)
...
subroutine, public create_atom_type(atom)
...
Definition atom_types.F:958
integer, parameter, public cgto_basis
Definition atom_types.F:69
integer, parameter, public gto_basis
Definition atom_types.F:69
integer, parameter, public sto_basis
Definition atom_types.F:69
subroutine, public release_atom_type(atom)
...
Definition atom_types.F:982
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.
Definition atom_types.F:378
logical function, public atom_compare_grids(grid1, grid2)
...
subroutine, public release_opgrid(opgrid)
...
subroutine, public release_atom_orbs(orbs)
...
integer, parameter, public lmat
Definition atom_types.F:67
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)
...
Definition atom_types.F:924
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
subroutine, public init_atom_basis_default_pp(basis)
...
Definition atom_types.F:715
subroutine, public create_opmat(opmat, n, lmax)
...
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
Definition atom_types.F:792
subroutine, public read_ecp_potential(element_symbol, potential, pseudo_name, pseudo_file, potential_section)
...
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.
Definition atom_upf.F:14
subroutine, public atom_read_upf(pot, upf_filename, read_header)
...
Definition atom_upf.F:102
pure subroutine, public atom_release_upf(upfpot)
...
Definition atom_upf.F:875
Definition atom.F:9
Calculates Bessel functions.
Definition bessel_lib.F:16
elemental impure real(kind=dp) function, public bessel0(x, l)
...
Definition bessel_lib.F:161
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public limpanuparb2011
logical function, public cp_sll_val_next(iterator, el_att)
returns true if the actual element is valid (i.e. iterator ont at end) moves the iterator to the next...
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.
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_rks_atom
integer, parameter, public do_analytic
integer, parameter, public sgp_pseudo
integer, parameter, public gth_pseudo
integer, parameter, public ecp_pseudo
integer, parameter, public do_nonrel_atom
integer, parameter, public no_conf
integer, parameter, public do_potential_mix_cl
integer, parameter, public upf_pseudo
integer, parameter, public contracted_gto
integer, parameter, public poly_conf
integer, parameter, public no_pseudo
integer, parameter, public barrier_conf
integer, parameter, public do_numeric
integer, parameter, public do_gapw_log
integer, parameter, public do_potential_coulomb
integer, parameter, public gaussian
integer, parameter, public do_potential_short
integer, parameter, public do_semi_analytic
integer, parameter, public geometrical_gto
integer, parameter, public do_potential_long
integer, parameter, public numerical
integer, parameter, public slater
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_list_get(section_vals, keyword_name, i_rep_section, list)
returns the requested list
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
a wrapper for basic fortran types.
subroutine, public val_get(val, has_l, has_i, has_r, has_lc, has_c, l_val, l_vals, i_val, i_vals, r_val, r_vals, c_val, c_vals, len_c, type_of_var, enum)
returns the stored values
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
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.
Definition atom_types.F:78
Provides all information about a pseudopotential.
Definition atom_types.F:98
Provides info about hartree-fock exchange (For now, we only support potentials that can be represente...
Definition atom_types.F:187
Information on optimization procedure.
Definition atom_types.F:283
Holds atomic orbitals and energies.
Definition atom_types.F:237
Provides all information on states and occupation.
Definition atom_types.F:198
Provides all information about an atomic kind.
Definition atom_types.F:293
Holds atomic integrals.
Definition atom_types.F:212
Operator grids.
Definition atom_types.F:258
Operator matrices.
Definition atom_types.F:251
represent a single linked list that stores pointers to the elements
a type to have a wrapper that stores any basic fortran type