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