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