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