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