(git:b77b4be)
Loading...
Searching...
No Matches
atom_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Define the atom type and its sub types
10!> \author jgh
11!> \date 03.03.2008
12!> \version 1.0
13!>
14! **************************************************************************************************
15MODULE atom_types
16 USE atom_upf, ONLY: atom_read_upf,&
17 atom_release_upf,&
18 atom_upfpot_type
19 USE bessel_lib, ONLY: bessel0
20 USE bibliography, ONLY: limpanuparb2011,&
21 cite_reference
22 USE cp_linked_list_input, ONLY: cp_sll_val_next,&
23 cp_sll_val_type
24 USE cp_parser_methods, ONLY: parser_get_next_line,&
25 parser_get_object,&
26 parser_read_line,&
27 parser_search_string,&
28 parser_test_next_token
29 USE cp_parser_types, ONLY: cp_parser_type,&
30 parser_create,&
31 parser_release
32 USE input_constants, ONLY: &
33 barrier_conf, contracted_gto, do_analytic, do_gapw_log, do_nonrel_atom, do_numeric, &
34 do_potential_coulomb, do_potential_long, do_potential_mix_cl, do_potential_short, &
35 do_rks_atom, do_semi_analytic, ecp_pseudo, gaussian, geometrical_gto, gth_pseudo, no_conf, &
36 no_pseudo, numerical, poly_conf, sgp_pseudo, slater, upf_pseudo
37 USE input_section_types, ONLY: section_vals_get,&
38 section_vals_get_subs_vals,&
39 section_vals_list_get,&
40 section_vals_type,&
41 section_vals_val_get
42 USE input_val_types, ONLY: val_get,&
43 val_type
44 USE kinds, ONLY: default_string_length,&
45 dp
46 USE mathconstants, ONLY: dfac,&
47 fac,&
48 pi,&
49 rootpi
50 USE periodic_table, ONLY: get_ptable_info,&
51 ptable
52 USE qs_grid_atom, ONLY: allocate_grid_atom,&
53 create_grid_atom,&
54 deallocate_grid_atom,&
55 grid_atom_type
56 USE string_utilities, ONLY: remove_word,&
57 uppercase
58#include "./base/base_uses.f90"
59
60 IMPLICIT NONE
61
62 PRIVATE
63
64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_types'
65
66 ! maximum l-quantum number considered in atomic code/basis
67 INTEGER, PARAMETER :: lmat = 5
68
69 INTEGER, PARAMETER :: gto_basis = 100, &
70 cgto_basis = 101, &
71 sto_basis = 102, &
72 num_basis = 103
73
74 INTEGER, PARAMETER :: nmax = 25
75
76!> \brief Provides all information about a basis set
77! **************************************************************************************************
78 TYPE atom_basis_type
79 INTEGER :: basis_type = gto_basis
80 INTEGER, DIMENSION(0:lmat) :: nbas = 0
81 INTEGER, DIMENSION(0:lmat) :: nprim = 0
82 REAL(kind=dp), DIMENSION(:, :), POINTER :: am => null() !GTO exponents
83 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: cm => null() !Contraction coeffs
84 REAL(kind=dp), DIMENSION(:, :), POINTER :: as => null() !STO exponents
85 INTEGER, DIMENSION(:, :), POINTER :: ns => null() !STO n-quantum numbers
86 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: bf => null() !num. bsf
87 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dbf => null() !derivatives (num)
88 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ddbf => null() !2nd derivatives (num)
89 REAL(kind=dp) :: eps_eig = 0.0_dp
90 TYPE(grid_atom_type), POINTER :: grid => null()
91 LOGICAL :: geometrical = .false.
92 REAL(kind=dp) :: aval = 0.0_dp, cval = 0.0_dp
93 INTEGER, DIMENSION(0:lmat) :: start = 0
94 END TYPE atom_basis_type
95
96!> \brief Provides all information about a pseudopotential
97! **************************************************************************************************
98 TYPE atom_gthpot_type
99 CHARACTER(LEN=2) :: symbol = ""
100 CHARACTER(LEN=default_string_length) :: pname = ""
101 INTEGER, DIMENSION(0:lmat) :: econf = 0
102 REAL(dp) :: zion = 0.0_dp
103 REAL(dp) :: rc = 0.0_dp
104 INTEGER :: ncl = 0
105 REAL(dp), DIMENSION(5) :: cl = 0.0_dp
106 INTEGER, DIMENSION(0:lmat) :: nl = 0
107 REAL(dp), DIMENSION(0:lmat) :: rcnl = 0.0_dp
108 REAL(dp), DIMENSION(4, 4, 0:lmat) :: hnl = 0.0_dp
109 ! 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
344CONTAINS
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 (0)
2984 CASE (2)
2985 potential%econf(0:3) = potential%econf(0:3) - ptable(2)%e_conv(0:3)
2986 CASE (10)
2987 potential%econf(0:3) = potential%econf(0:3) - ptable(10)%e_conv(0:3)
2988 CASE (18)
2989 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
2990 CASE (28)
2991 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
2992 potential%econf(2) = potential%econf(2) - 10
2993 CASE (36)
2994 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
2995 CASE (46)
2996 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
2997 potential%econf(2) = potential%econf(2) - 10
2998 CASE (54)
2999 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3000 CASE (60)
3001 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3002 potential%econf(2) = potential%econf(2) - 10
3003 potential%econf(3) = potential%econf(3) - 14
3004 CASE (68)
3005 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3006 potential%econf(3) = potential%econf(3) - 14
3007 CASE (78)
3008 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3009 potential%econf(2) = potential%econf(2) - 10
3010 potential%econf(3) = potential%econf(3) - 14
3011 END SELECT
3012 !
3013 cpassert(all(potential%econf >= 0))
3014
3015 END SUBROUTINE read_ecp_potential
3016! **************************************************************************************************
3017!> \brief ...
3018!> \param grid1 ...
3019!> \param grid2 ...
3020!> \return ...
3021! **************************************************************************************************
3022 FUNCTION atom_compare_grids(grid1, grid2) RESULT(is_equal)
3023 TYPE(grid_atom_type) :: grid1, grid2
3024 LOGICAL :: is_equal
3025
3026 INTEGER :: i
3027 REAL(kind=dp) :: dr, dw
3028
3029 is_equal = .true.
3030 IF (grid1%nr == grid2%nr) THEN
3031 DO i = 1, grid2%nr
3032 dr = abs(grid1%rad(i) - grid2%rad(i))
3033 dw = abs(grid1%wr(i) - grid2%wr(i))
3034 IF (dr + dw > 1.0e-12_dp) THEN
3035 is_equal = .false.
3036 EXIT
3037 END IF
3038 END DO
3039 ELSE
3040 is_equal = .false.
3041 END IF
3042
3043 END FUNCTION atom_compare_grids
3044! **************************************************************************************************
3045
3046END MODULE atom_types
cp_parser_methods::parser_get_object
Definition cp_parser_methods.F:54
atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
atom_types::num_basis
integer, parameter, public num_basis
Definition atom_types.F:69
atom_types::read_atom_opt_section
subroutine, public read_atom_opt_section(optimization, opt_section)
...
Definition atom_types.F: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:3023
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:21
bibliography::limpanuparb2011
integer, save, public limpanuparb2011
Definition bibliography.F:36
cp_linked_list_input
Definition cp_linked_list_input.F:9
cp_linked_list_input::cp_sll_val_next
logical function, public cp_sll_val_next(iterator, el_att)
returns true if the actual element is valid (i.e. iterator ont at end) moves the iterator to the next...
Definition cp_linked_list_input.F:1474
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::parser_read_line
subroutine, public parser_read_line(parser, nline, at_end)
Read the next line from a logical unit "unit" (I/O node only). Skip (nline-1) lines and skip also all...
Definition cp_parser_methods.F:147
cp_parser_methods::parser_get_next_line
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Definition cp_parser_methods.F:421
cp_parser_methods::parser_test_next_token
character(len=3) function, public parser_test_next_token(parser, string_length)
Test next input object.
Definition cp_parser_methods.F:691
cp_parser_methods::parser_search_string
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Definition cp_parser_methods.F:786
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition cp_parser_types.F:109
cp_parser_types::parser_create
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Definition cp_parser_types.F:147
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_rks_atom
integer, parameter, public do_rks_atom
Definition input_constants.F:159
input_constants::do_analytic
integer, parameter, public do_analytic
Definition input_constants.F:149
input_constants::sgp_pseudo
integer, parameter, public sgp_pseudo
Definition input_constants.F:174
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:174
input_constants::ecp_pseudo
integer, parameter, public ecp_pseudo
Definition input_constants.F:174
input_constants::do_nonrel_atom
integer, parameter, public do_nonrel_atom
Definition input_constants.F:165
input_constants::no_conf
integer, parameter, public no_conf
Definition input_constants.F:180
input_constants::do_potential_mix_cl
integer, parameter, public do_potential_mix_cl
Definition input_constants.F:807
input_constants::upf_pseudo
integer, parameter, public upf_pseudo
Definition input_constants.F:174
input_constants::contracted_gto
integer, parameter, public contracted_gto
Definition input_constants.F:153
input_constants::poly_conf
integer, parameter, public poly_conf
Definition input_constants.F:180
input_constants::no_pseudo
integer, parameter, public no_pseudo
Definition input_constants.F:174
input_constants::barrier_conf
integer, parameter, public barrier_conf
Definition input_constants.F:180
input_constants::do_numeric
integer, parameter, public do_numeric
Definition input_constants.F:149
input_constants::do_gapw_log
integer, parameter, public do_gapw_log
Definition input_constants.F:226
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:807
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:807
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:807
input_constants::numerical
integer, parameter, public numerical
Definition input_constants.F:153
input_constants::slater
integer, parameter, public slater
Definition input_constants.F:153
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_list_get
subroutine, public section_vals_list_get(section_vals, keyword_name, i_rep_section, list)
returns the requested list
Definition input_section_types.F:1156
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
input_val_types
a wrapper for basic fortran types.
Definition input_val_types.F:14
input_val_types::val_get
subroutine, public val_get(val, has_l, has_i, has_r, has_lc, has_c, l_val, l_vals, i_val, i_vals, r_val, r_vals, c_val, c_vals, len_c, type_of_var, enum)
returns the stored values
Definition input_val_types.F:334
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
qs_grid_atom
Definition qs_grid_atom.F:8
qs_grid_atom::deallocate_grid_atom
subroutine, public deallocate_grid_atom(grid_atom)
Deallocate a Gaussian-type orbital (GTO) basis set data set.
Definition qs_grid_atom.F:105
qs_grid_atom::allocate_grid_atom
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
Definition qs_grid_atom.F:73
qs_grid_atom::create_grid_atom
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Definition qs_grid_atom.F:183
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::newline
character(len=1), parameter, public newline
Definition string_utilities.F:29
string_utilities::remove_word
subroutine, public remove_word(string)
remove a word from a string (words are separated by white spaces)
Definition string_utilities.F:3275
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atom_types::atom_basis_type
Provides all information about a basis set.
Definition atom_types.F:78
atom_types::atom_ecppot_type
Definition atom_types.F:130
atom_types::atom_gthpot_type
Provides all information about a pseudopotential.
Definition atom_types.F:98
atom_types::atom_hfx_type
Provides info about hartree-fock exchange (For now, we only support potentials that can be represente...
Definition atom_types.F:184
atom_types::atom_integrals
Definition atom_types.F:213
atom_types::atom_optimization_type
Information on optimization procedure.
Definition atom_types.F:280
atom_types::atom_orbitals
Holds atomic orbitals and energies.
Definition atom_types.F:234
atom_types::atom_p_type
Definition atom_types.F:320
atom_types::atom_potential_type
Definition atom_types.F:168
atom_types::atom_sgppot_type
Definition atom_types.F:146
atom_types::atom_state
Provides all information on states and occupation.
Definition atom_types.F:195
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:290
atom_types::eri
Holds atomic integrals.
Definition atom_types.F:209
atom_types::opgrid_type
Operator grids.
Definition atom_types.F:255
atom_types::opmat_type
Operator matrices.
Definition atom_types.F:248
atom_upf::atom_upfpot_type
Definition atom_upf.F:32
cp_linked_list_input::cp_sll_val_type
represent a single linked list that stores pointers to the elements
Definition cp_linked_list_input.F:341
cp_parser_types::cp_parser_type
represent a parser
Definition cp_parser_types.F:77
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
input_val_types::val_type
a type to have a wrapper that stores any basic fortran type
Definition input_val_types.F:63
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44