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