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