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