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 Read xTB parameters.
10 : !> \author JGH (10.2018)
11 : ! **************************************************************************************************
12 : MODULE xtb_parameters
13 :
14 : USE basis_set_types, ONLY: allocate_sto_basis_set,&
15 : create_gto_from_sto_basis,&
16 : deallocate_sto_basis_set,&
17 : gto_basis_set_type,&
18 : set_sto_basis_set,&
19 : sto_basis_set_type
20 : USE cp_control_types, ONLY: xtb_control_type
21 : USE cp_linked_list_input, ONLY: cp_sll_val_next,&
22 : cp_sll_val_type
23 : USE cp_parser_methods, ONLY: parser_get_next_line,&
24 : parser_get_object
25 : USE cp_parser_types, ONLY: cp_parser_type,&
26 : parser_create,&
27 : parser_release
28 : USE input_section_types, ONLY: section_vals_get,&
29 : section_vals_get_subs_vals,&
30 : section_vals_list_get,&
31 : section_vals_type
32 : USE input_val_types, ONLY: val_get,&
33 : val_type
34 : USE kinds, ONLY: default_string_length,&
35 : dp
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE periodic_table, ONLY: get_ptable_info,&
38 : ptable
39 : USE physcon, ONLY: bohr,&
40 : evolt
41 : USE string_utilities, ONLY: remove_word,&
42 : uppercase
43 : USE xtb_types, ONLY: xtb_atom_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 :
48 : PRIVATE
49 :
50 : INTEGER, PARAMETER, PRIVATE :: nelem = 106
51 : ! H He
52 : ! Li Be B C N O F Ne
53 : ! Na Mg Al Si P S Cl Ar
54 : ! K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
55 : ! Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe
56 : ! Cs Ba La Ce-Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
57 : ! Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Ha 106
58 :
59 : !&<
60 : ! Element Valence
61 : INTEGER, DIMENSION(0:nelem), &
62 : PARAMETER, PRIVATE :: zval = (/-1, & ! 0
63 : 1, 2, & ! 2
64 : 1, 2, 3, 4, 5, 6, 7, 8, & ! 10
65 : 1, 2, 3, 4, 5, 6, 7, 8, & ! 18
66 : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 36
67 : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 54
68 : 1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, &
69 : 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 86
70 : -1, -1, -1, 4, -1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
71 : !&>
72 :
73 : !&<
74 : ! Element Pauling Electronegativity
75 : REAL(KIND=dp), DIMENSION(0:nelem), &
76 : PARAMETER, PRIVATE :: eneg = (/0.00_dp, & ! 0
77 : 2.20_dp, 3.00_dp, & ! 2
78 : 0.98_dp, 1.57_dp, 2.04_dp, 2.55_dp, 3.04_dp, 3.44_dp, 3.98_dp, 4.50_dp, & ! 10
79 : 0.93_dp, 1.31_dp, 1.61_dp, 1.90_dp, 2.19_dp, 2.58_dp, 3.16_dp, 3.50_dp, & ! 18
80 : 0.82_dp, 1.00_dp, 1.36_dp, 1.54_dp, 1.63_dp, 1.66_dp, 1.55_dp, 1.83_dp, &
81 : 1.88_dp, 1.91_dp, 1.90_dp, 1.65_dp, 1.81_dp, 2.01_dp, 2.18_dp, 2.55_dp, 2.96_dp, 3.00_dp, & ! 36
82 : 0.82_dp, 0.95_dp, 1.22_dp, 1.33_dp, 1.60_dp, 2.16_dp, 1.90_dp, 2.20_dp, &
83 : 2.28_dp, 2.20_dp, 1.93_dp, 1.69_dp, 1.78_dp, 1.96_dp, 2.05_dp, 2.10_dp, 2.66_dp, 2.60_dp, & ! 54
84 : 0.79_dp, 0.89_dp, 1.10_dp, &
85 : 1.12_dp, 1.13_dp, 1.14_dp, 1.15_dp, 1.17_dp, 1.18_dp, 1.20_dp, 1.21_dp, &
86 : 1.22_dp, 1.23_dp, 1.24_dp, 1.25_dp, 1.26_dp, 1.27_dp, & ! Lanthanides
87 : 1.30_dp, 1.50_dp, 2.36_dp, 1.90_dp, 2.20_dp, 2.20_dp, 2.28_dp, 2.54_dp, &
88 : 2.00_dp, 2.04_dp, 2.33_dp, 2.02_dp, 2.00_dp, 2.20_dp, 2.20_dp, & ! 86
89 : 0.70_dp, 0.89_dp, 1.10_dp, &
90 : 1.30_dp, 1.50_dp, 1.38_dp, 1.36_dp, 1.28_dp, 1.30_dp, 1.30_dp, 1.30_dp, &
91 : 1.30_dp, 1.30_dp, 1.30_dp, 1.30_dp, 1.30_dp, 1.50_dp, & ! Actinides
92 : 1.50_dp, 1.50_dp, 1.50_dp/)
93 : !&>
94 :
95 : !&<
96 : ! Shell occupation
97 : INTEGER, DIMENSION(1:5, 0:nelem) :: occupation = RESHAPE((/0,0,0,0,0, & ! 0
98 : 1,0,0,0,0, 2,0,0,0,0, & ! 2
99 : 1,0,0,0,0, 2,0,0,0,0, 2,1,0,0,0, 2,2,0,0,0, 2,3,0,0,0, 2,4,0,0,0, 2,5,0,0,0, 2,6,0,0,0, & ! 10
100 : 1,0,0,0,0, 2,0,0,0,0, 2,1,0,0,0, 2,2,0,0,0, 2,3,0,0,0, 2,4,0,0,0, 2,5,0,0,0, 2,6,0,0,0, & ! 18
101 : 1,0,0,0,0, 2,0,0,0,0, 2,0,1,0,0, 2,0,2,0,0, 2,0,3,0,0, 2,0,4,0,0, 2,0,5,0,0, 2,0,6,0,0, &
102 : 2,0,7,0,0, 2,0,8,0,0, 2,0,9,0,0, 2,0,0,0,0, 2,1,0,0,0, 2,2,0,0,0, 2,3,0,0,0, 2,4,0,0,0, 2,5,0,0,0, 2,6,0,0,0, & ! 36
103 : 1,0,0,0,0, 2,0,0,0,0, 2,0,1,0,0, 2,0,2,0,0, 2,0,3,0,0, 2,0,4,0,0, 2,0,5,0,0, 2,0,6,0,0, & !
104 : 2,0,7,0,0, 2,0,8,0,0, 2,0,9,0,0, 2,0,0,0,0, 2,1,0,0,0, 2,2,0,0,0, 2,3,0,0,0, 2,4,0,0,0, 2,5,0,0,0, 2,6,0,0,0, & ! 54
105 : 1,0,0,0,0, 2,0,0,0,0, 2,0,1,0,0, &
106 : 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, &
107 : 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, 2,0,1,0,0, & ! Lanthanides
108 : 2,0,2,0,0, 2,0,3,0,0, 2,0,4,0,0, 2,0,5,0,0, 2,0,6,0,0, 2,0,7,0,0, 2,0,8,0,0, 2,0,9,0,0, &
109 : 2,0,0,0,0, 2,1,0,0,0, 2,2,0,0,0, 2,3,0,0,0, 2,4,0,0,0, 2,5,0,0,0, 2,6,0,0,0, & ! 86 (last element defined)
110 : 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, & !
111 : 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, &
112 : 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, & ! Actinides
113 : 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0/), (/5, nelem+1/))
114 : !&>
115 :
116 : !&<
117 : ! COVALENT RADII
118 : ! based on "Atomic Radii of the Elements," M. Mantina, R. Valero, C. J. Cramer, and D. G. Truhlar,
119 : ! in CRC Handbook of Chemistry and Physics, 91st Edition (2010-2011),
120 : ! edited by W. M. Haynes (CRC Press, Boca Raton, FL, 2010), pages 9-49-9-50;
121 : ! corrected Nov. 17, 2010 for the 92nd edition.
122 : REAL(KIND=dp), DIMENSION(0:nelem), &
123 : PARAMETER, PRIVATE :: crad = (/0.00_dp, & ! 0
124 : 0.32_dp, 0.37_dp, & ! 2
125 : 1.30_dp, 0.99_dp, 0.84_dp, 0.75_dp, 0.71_dp, 0.64_dp, 0.60_dp, 0.62_dp, & ! 10
126 : 1.60_dp, 1.40_dp, 1.24_dp, 1.14_dp, 1.09_dp, 1.04_dp, 1.00_dp, 1.01_dp, & ! 18
127 : 2.00_dp, 1.74_dp, 1.59_dp, 1.48_dp, 1.44_dp, 1.30_dp, 1.29_dp, 1.24_dp, &
128 : 1.18_dp, 1.17_dp, 1.22_dp, 1.20_dp, 1.23_dp, 1.20_dp, 1.20_dp, 1.18_dp, 1.17_dp, 1.16_dp, & ! 36
129 : 2.15_dp, 1.90_dp, 1.76_dp, 1.64_dp, 1.56_dp, 1.46_dp, 1.38_dp, 1.36_dp, &
130 : 1.34_dp, 1.30_dp, 1.36_dp, 1.40_dp, 1.42_dp, 1.40_dp, 1.40_dp, 1.37_dp, 1.36_dp, 1.36_dp, & ! 54
131 : 2.38_dp, 2.06_dp, 1.94_dp, &
132 : 1.84_dp, 1.90_dp, 1.88_dp, 1.86_dp, 1.85_dp, 1.83_dp, 1.82_dp, 1.81_dp, &
133 : 1.80_dp, 1.79_dp, 1.77_dp, 1.77_dp, 1.78_dp, 1.74_dp, & ! Lanthanides
134 : 1.64_dp, 1.58_dp, 1.50_dp, 1.41_dp, 1.36_dp, 1.32_dp, 1.30_dp, 1.30_dp, &
135 : 1.32_dp, 1.44_dp, 1.45_dp, 1.50_dp, 1.42_dp, 1.48_dp, 1.46_dp, & ! 86
136 : 2.42_dp, 2.11_dp, 2.01_dp, &
137 : 1.90_dp, 1.84_dp, 1.83_dp, 1.80_dp, 1.80_dp, 1.51_dp, 0.96_dp, 1.54_dp, &
138 : 1.83_dp, 1.50_dp, 1.50_dp, 1.50_dp, 1.50_dp, 1.50_dp, & ! Actinides
139 : 1.50_dp, 1.50_dp, 1.50_dp/)
140 : !&>
141 :
142 : !&<
143 : ! Charge Limits (Mulliken)
144 : REAL(KIND=dp), DIMENSION(0:nelem), &
145 : PARAMETER, PRIVATE :: clmt = (/0.00_dp, & ! 0
146 : 1.05_dp, 1.25_dp, & ! 2
147 : 1.05_dp, 2.05_dp, 3.00_dp, 4.00_dp, 3.00_dp, 2.00_dp, 1.25_dp, 1.00_dp, & ! 10
148 : 1.05_dp, 2.05_dp, 3.00_dp, 4.00_dp, 3.00_dp, 2.00_dp, 1.25_dp, 1.00_dp, & ! 18
149 : 1.05_dp, 2.05_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
150 : 3.50_dp, 3.50_dp, 3.50_dp, 2.50_dp, 2.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 1.25_dp, 1.00_dp, & ! 36
151 : 1.05_dp, 2.05_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
152 : 3.50_dp, 3.50_dp, 3.50_dp, 2.50_dp, 2.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 1.25_dp, 1.00_dp, & ! 54
153 : 1.05_dp, 2.05_dp, 3.00_dp, &
154 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, &
155 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, & ! Lanthanides
156 : 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
157 : 2.50_dp, 2.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 1.25_dp, 1.00_dp, & ! 86
158 : 1.05_dp, 2.05_dp, 3.00_dp, &
159 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, &
160 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, & ! Actinides
161 : 3.00_dp, 3.00_dp, 3.00_dp/)
162 : !&>
163 :
164 : ! *** Global parameters ***
165 :
166 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_parameters'
167 :
168 : ! *** Public data types ***
169 :
170 : PUBLIC :: xtb_parameters_init, xtb_parameters_read, xtb_parameters_set, init_xtb_basis, &
171 : xtb_set_kab
172 :
173 : CONTAINS
174 :
175 : ! **************************************************************************************************
176 : !> \brief ...
177 : !> \param param ...
178 : !> \param element_symbol ...
179 : !> \param parameter_file_path ...
180 : !> \param parameter_file_name ...
181 : !> \param para_env ...
182 : ! **************************************************************************************************
183 1184 : SUBROUTINE xtb_parameters_init(param, element_symbol, parameter_file_path, parameter_file_name, &
184 : para_env)
185 :
186 : TYPE(xtb_atom_type), POINTER :: param
187 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
188 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
189 : TYPE(mp_para_env_type), POINTER :: para_env
190 :
191 : CHARACTER(len=2) :: enam, esym
192 : CHARACTER(len=default_string_length) :: aname, filename
193 : INTEGER :: i, ia, l
194 : LOGICAL :: at_end, found
195 : TYPE(cp_parser_type) :: parser
196 :
197 592 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(parameter_file_name))
198 592 : CALL parser_create(parser, filename, para_env=para_env)
199 592 : found = .FALSE.
200 : DO
201 : at_end = .FALSE.
202 9856 : CALL parser_get_next_line(parser, 1, at_end)
203 9856 : IF (at_end) EXIT
204 9854 : CALL parser_get_object(parser, aname)
205 9854 : enam = aname
206 9854 : esym = element_symbol
207 9854 : CALL uppercase(enam)
208 9854 : CALL uppercase(esym)
209 9856 : IF (enam == esym) THEN
210 590 : found = .TRUE.
211 590 : CALL parser_get_object(parser, param%eta)
212 590 : CALL parser_get_object(parser, param%xgamma)
213 590 : CALL parser_get_object(parser, param%alpha)
214 590 : CALL parser_get_object(parser, param%zneff)
215 1844 : DO i = 1, 5
216 1844 : CALL parser_get_object(parser, aname)
217 1844 : ia = ICHAR(aname(1:1))
218 1844 : IF (ia >= 49 .AND. ia <= 57) THEN
219 1254 : CALL parser_get_object(parser, param%kpoly(i))
220 1254 : CALL parser_get_object(parser, param%kappa(i))
221 1254 : CALL parser_get_object(parser, param%hen(i))
222 1254 : CALL parser_get_object(parser, param%zeta(i))
223 1254 : param%nshell = i
224 1254 : param%nval(i) = ia - 48
225 816 : SELECT CASE (aname(2:2))
226 : CASE ("s", "S")
227 816 : param%lval(i) = 0
228 : CASE ("p", "P")
229 364 : param%lval(i) = 1
230 : CASE ("d", "D")
231 74 : param%lval(i) = 2
232 : CASE ("f", "F")
233 0 : param%lval(i) = 3
234 : CASE DEFAULT
235 1254 : CPABORT("xTB PARAMETER ERROR")
236 : END SELECT
237 1254 : CALL parser_get_next_line(parser, 1, at_end)
238 1254 : IF (at_end) EXIT
239 : ELSE
240 : EXIT
241 : END IF
242 : END DO
243 : EXIT
244 : END IF
245 : END DO
246 592 : IF (found) THEN
247 590 : param%typ = "STANDARD"
248 590 : param%symbol = element_symbol
249 590 : param%defined = .TRUE.
250 590 : CALL get_ptable_info(element_symbol, number=ia)
251 590 : param%z = ia
252 590 : param%aname = ptable(ia)%name
253 1844 : param%lmax = MAXVAL(param%lval(1:param%nshell))
254 590 : param%natorb = 0
255 1844 : DO i = 1, param%nshell
256 1254 : l = param%lval(i)
257 1844 : param%natorb = param%natorb + (2*l + 1)
258 : END DO
259 590 : param%zeff = zval(ia)
260 : ELSE
261 2 : esym = element_symbol
262 2 : CALL uppercase(esym)
263 2 : IF ("X " == esym) THEN
264 2 : param%typ = "GHOST"
265 2 : param%symbol = element_symbol
266 2 : param%defined = .FALSE.
267 2 : param%z = 0
268 2 : param%aname = "X "
269 2 : param%lmax = 0
270 2 : param%natorb = 0
271 2 : param%nshell = 0
272 2 : param%zeff = 0.0_dp
273 : ELSE
274 0 : param%defined = .FALSE.
275 : CALL cp_warn(__LOCATION__, "xTB parameters for element "//element_symbol// &
276 0 : " were not found in the parameter file "//ADJUSTL(TRIM(filename)))
277 : END IF
278 : END IF
279 592 : CALL parser_release(parser)
280 :
281 1776 : END SUBROUTINE xtb_parameters_init
282 :
283 : ! **************************************************************************************************
284 : !> \brief Read atom parameters for xTB Hamiltonian from input file
285 : !> \param param ...
286 : !> \param element_symbol ...
287 : !> \param xtb_section ...
288 : ! **************************************************************************************************
289 1184 : SUBROUTINE xtb_parameters_read(param, element_symbol, xtb_section)
290 :
291 : TYPE(xtb_atom_type), POINTER :: param
292 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
293 : TYPE(section_vals_type), POINTER :: xtb_section
294 :
295 : CHARACTER(LEN=2) :: label
296 : CHARACTER(len=20*default_string_length) :: line_att
297 : INTEGER :: i, ia, k, l, nshell
298 : LOGICAL :: explicit, found, is_ok
299 : TYPE(cp_sll_val_type), POINTER :: list
300 : TYPE(section_vals_type), POINTER :: ap_section
301 : TYPE(val_type), POINTER :: val
302 :
303 : !
304 : ! This could probably be done nicer
305 : !
306 592 : NULLIFY (list, val)
307 592 : ap_section => section_vals_get_subs_vals(xtb_section, "ATOM_PARAMETER")
308 592 : CALL section_vals_get(ap_section, explicit=explicit)
309 592 : IF (explicit) THEN
310 12 : CALL section_vals_list_get(ap_section, "_DEFAULT_KEYWORD_", list=list)
311 12 : found = .FALSE.
312 12 : nshell = 0
313 : DO
314 76 : is_ok = cp_sll_val_next(list, val)
315 76 : IF (.NOT. is_ok) EXIT
316 68 : CALL val_get(val, c_val=line_att)
317 80 : IF (found) THEN
318 12 : READ (line_att, *) label
319 12 : CALL remove_word(line_att)
320 12 : ia = ICHAR(label(1:1))
321 12 : IF (ia >= 49 .AND. ia <= 57) THEN
322 8 : nshell = nshell + 1
323 8 : k = nshell
324 8 : param%nval(k) = ia - 48
325 0 : SELECT CASE (label(2:2))
326 : CASE ("s", "S")
327 0 : param%lval(k) = 0
328 : CASE ("p", "P")
329 8 : param%lval(k) = 1
330 : CASE ("d", "D")
331 0 : param%lval(k) = 2
332 : CASE ("f", "F")
333 0 : param%lval(k) = 3
334 : CASE DEFAULT
335 8 : CPABORT("xTB PARAMETER ERROR")
336 : END SELECT
337 : !
338 8 : READ (line_att, *) param%kpoly(k)
339 8 : CALL remove_word(line_att)
340 8 : READ (line_att, *) param%kappa(k)
341 8 : CALL remove_word(line_att)
342 8 : READ (line_att, *) param%hen(k)
343 8 : CALL remove_word(line_att)
344 8 : READ (line_att, *) param%zeta(k)
345 8 : CALL remove_word(line_att)
346 : ELSE
347 : EXIT
348 : END IF
349 : ELSE
350 56 : READ (line_att, *) label
351 56 : CALL remove_word(line_att)
352 56 : IF (label == element_symbol) THEN
353 8 : found = .TRUE.
354 8 : nshell = nshell + 1
355 8 : k = nshell
356 8 : READ (line_att, *) param%eta
357 8 : CALL remove_word(line_att)
358 8 : READ (line_att, *) param%xgamma
359 8 : CALL remove_word(line_att)
360 8 : READ (line_att, *) param%alpha
361 8 : CALL remove_word(line_att)
362 8 : READ (line_att, *) param%zneff
363 8 : CALL remove_word(line_att)
364 8 : READ (line_att, *) label
365 8 : CALL remove_word(line_att)
366 8 : ia = ICHAR(label(1:1))
367 8 : CPASSERT((ia >= 49 .AND. ia <= 57))
368 8 : param%nval(k) = ia - 48
369 8 : SELECT CASE (label(2:2))
370 : CASE ("s", "S")
371 8 : param%lval(k) = 0
372 : CASE ("p", "P")
373 0 : param%lval(k) = 1
374 : CASE ("d", "D")
375 0 : param%lval(k) = 2
376 : CASE ("f", "F")
377 0 : param%lval(k) = 3
378 : CASE DEFAULT
379 8 : CPABORT("xTB PARAMETER ERROR")
380 : END SELECT
381 : !
382 8 : READ (line_att, *) param%kpoly(k)
383 8 : CALL remove_word(line_att)
384 8 : READ (line_att, *) param%kappa(k)
385 8 : CALL remove_word(line_att)
386 8 : READ (line_att, *) param%hen(k)
387 8 : CALL remove_word(line_att)
388 8 : READ (line_att, *) param%zeta(k)
389 8 : CALL remove_word(line_att)
390 : END IF
391 : END IF
392 : END DO
393 12 : IF (found) THEN
394 8 : param%typ = "STANDARD"
395 8 : param%symbol = element_symbol
396 8 : param%defined = .TRUE.
397 8 : CALL get_ptable_info(element_symbol, number=ia)
398 8 : param%z = ia
399 8 : param%aname = ptable(ia)%name
400 24 : param%lmax = MAXVAL(param%lval(1:param%nshell))
401 8 : param%natorb = 0
402 8 : param%nshell = nshell
403 24 : DO i = 1, param%nshell
404 16 : l = param%lval(i)
405 24 : param%natorb = param%natorb + (2*l + 1)
406 : END DO
407 8 : param%zeff = zval(ia)
408 : END IF
409 : END IF
410 :
411 592 : END SUBROUTINE xtb_parameters_read
412 :
413 : ! **************************************************************************************************
414 : !> \brief Read atom parameters for xTB Hamiltonian from input file
415 : !> \param param ...
416 : ! **************************************************************************************************
417 592 : SUBROUTINE xtb_parameters_set(param)
418 :
419 : TYPE(xtb_atom_type), POINTER :: param
420 :
421 : INTEGER :: i, is, l, na
422 : REAL(KIND=dp), DIMENSION(5) :: kp
423 :
424 592 : IF (param%defined) THEN
425 : ! AO to shell pointer
426 : ! AO to l-qn pointer
427 590 : na = 0
428 1844 : DO is = 1, param%nshell
429 1254 : l = param%lval(is)
430 4122 : DO i = 1, 2*l + 1
431 2278 : na = na + 1
432 2278 : param%nao(na) = is
433 3532 : param%lao(na) = l
434 : END DO
435 : END DO
436 : !
437 590 : i = param%z
438 : ! Electronegativity
439 590 : param%electronegativity = eneg(i)
440 : ! covalent radius
441 590 : param%rcov = crad(i)*bohr
442 : ! shell occupances
443 3540 : param%occupation(:) = occupation(:, i)
444 : ! check for consistency
445 3540 : IF (ABS(param%zeff - SUM(param%occupation)) > 1.E-10_dp) THEN
446 0 : CALL cp_abort(__LOCATION__, "Element <"//TRIM(param%aname)//"> has inconsistent shell occupations")
447 : END IF
448 : ! orbital energies [evolt] -> [a.u.]
449 3540 : param%hen = param%hen/evolt
450 : ! some forgotten scaling parameters (not in orig. paper)
451 590 : param%xgamma = 0.1_dp*param%xgamma
452 3540 : param%kpoly(:) = 0.01_dp*param%kpoly(:)
453 3540 : param%kappa(:) = 0.1_dp*param%kappa(:)
454 : ! we have 1/6 g * q**3 (not 1/3)
455 590 : param%xgamma = -2.0_dp*param%xgamma
456 : ! we need kappa l-indexed
457 3540 : kp(:) = param%kappa(:)
458 3540 : param%kappa(:) = 0.0_dp
459 1844 : DO is = 1, param%nshell
460 1254 : l = param%lval(is)
461 1844 : IF (param%kappa(l + 1) == 0.0_dp) THEN
462 1254 : param%kappa(l + 1) = kp(is)
463 : ELSE
464 0 : CPASSERT(ABS(param%kappa(l + 1) - kp(is)) < 1.e-10_dp)
465 : END IF
466 : END DO
467 : ! kx
468 590 : IF (param%kx < -10._dp) THEN
469 : ! use defaults
470 1168 : SELECT CASE (param%z)
471 : CASE DEFAULT
472 578 : param%kx = 0.0_dp
473 : CASE (35) ! Br
474 12 : param%kx = 0.1_dp*0.381742_dp
475 : CASE (53) ! I
476 0 : param%kx = 0.1_dp*0.321944_dp
477 : CASE (85) ! At
478 590 : param%kx = 0.1_dp*0.220000_dp
479 : END SELECT
480 : END IF
481 : ! chmax
482 590 : param%chmax = clmt(i)
483 : END IF
484 :
485 592 : END SUBROUTINE xtb_parameters_set
486 :
487 : ! **************************************************************************************************
488 : !> \brief ...
489 : !> \param param ...
490 : !> \param gto_basis_set ...
491 : !> \param ngauss ...
492 : ! **************************************************************************************************
493 590 : SUBROUTINE init_xtb_basis(param, gto_basis_set, ngauss)
494 :
495 : TYPE(xtb_atom_type), POINTER :: param
496 : TYPE(gto_basis_set_type), POINTER :: gto_basis_set
497 : INTEGER, INTENT(IN) :: ngauss
498 :
499 590 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: symbol
500 : INTEGER :: i, nshell
501 590 : INTEGER, DIMENSION(:), POINTER :: lq, nq
502 590 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet
503 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
504 :
505 590 : IF (ASSOCIATED(param)) THEN
506 590 : IF (param%defined) THEN
507 590 : NULLIFY (sto_basis_set)
508 590 : CALL allocate_sto_basis_set(sto_basis_set)
509 590 : nshell = param%nshell
510 :
511 1770 : ALLOCATE (symbol(1:nshell))
512 1844 : symbol = ""
513 1844 : DO i = 1, nshell
514 590 : SELECT CASE (param%lval(i))
515 : CASE (0)
516 816 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "S"
517 : CASE (1)
518 364 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "P"
519 : CASE (2)
520 74 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "D"
521 : CASE (3)
522 0 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "F"
523 : CASE DEFAULT
524 1254 : CPABORT('BASIS SET OUT OF RANGE (lval)')
525 : END SELECT
526 : END DO
527 :
528 590 : IF (nshell > 0) THEN
529 3540 : ALLOCATE (nq(nshell), lq(nshell), zet(nshell))
530 3688 : nq(1:nshell) = param%nval(1:nshell)
531 3688 : lq(1:nshell) = param%lval(1:nshell)
532 3688 : zet(1:nshell) = param%zeta(1:nshell)
533 : CALL set_sto_basis_set(sto_basis_set, name=param%aname, nshell=nshell, symbol=symbol, &
534 590 : nq=nq, lq=lq, zet=zet)
535 590 : CALL create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss=ngauss, ortho=.TRUE.)
536 : END IF
537 :
538 : ! this will remove the allocated arrays
539 590 : CALL deallocate_sto_basis_set(sto_basis_set)
540 590 : DEALLOCATE (symbol, nq, lq, zet)
541 : END IF
542 :
543 : ELSE
544 0 : CPABORT("The pointer param is not associated")
545 : END IF
546 :
547 590 : END SUBROUTINE init_xtb_basis
548 :
549 : ! **************************************************************************************************
550 : !> \brief ...
551 : !> \param za ...
552 : !> \param zb ...
553 : !> \param xtb_control ...
554 : !> \return ...
555 : ! **************************************************************************************************
556 21816 : FUNCTION xtb_set_kab(za, zb, xtb_control) RESULT(kab)
557 :
558 : INTEGER, INTENT(IN) :: za, zb
559 : TYPE(xtb_control_type), INTENT(IN), POINTER :: xtb_control
560 : REAL(KIND=dp) :: kab
561 :
562 : INTEGER :: j, z
563 : LOGICAL :: custom
564 :
565 21816 : kab = 1.0_dp
566 21816 : custom = .FALSE.
567 :
568 21816 : IF (xtb_control%kab_nval .GT. 0) THEN
569 28 : DO j = 1, xtb_control%kab_nval
570 : IF ((za == xtb_control%kab_types(1, j) .AND. &
571 28 : zb == xtb_control%kab_types(2, j)) .OR. &
572 : (za == xtb_control%kab_types(2, j) .AND. &
573 0 : zb == xtb_control%kab_types(1, j))) THEN
574 28 : custom = .TRUE.
575 28 : kab = xtb_control%kab_vals(j)
576 28 : EXIT
577 : END IF
578 : END DO
579 : END IF
580 :
581 21816 : IF (.NOT. custom) THEN
582 21788 : IF (za == 1 .OR. zb == 1) THEN
583 : ! hydrogen
584 12472 : z = za + zb - 1
585 3100 : SELECT CASE (z)
586 : CASE (1)
587 3100 : kab = 0.96_dp
588 : CASE (5)
589 0 : kab = 0.95_dp
590 : CASE (7)
591 760 : kab = 1.04_dp
592 : CASE (28)
593 0 : kab = 0.90_dp
594 : CASE (75)
595 0 : kab = 0.80_dp
596 : CASE (78)
597 12472 : kab = 0.80_dp
598 : END SELECT
599 9316 : ELSEIF (za == 5 .OR. zb == 5) THEN
600 : ! Boron
601 0 : z = za + zb - 5
602 0 : SELECT CASE (z)
603 : CASE (15)
604 0 : kab = 0.97_dp
605 : END SELECT
606 9316 : ELSEIF (za == 7 .OR. zb == 7) THEN
607 : ! Nitrogen
608 2008 : z = za + zb - 7
609 0 : SELECT CASE (z)
610 : CASE (14)
611 : !xtb orig code parameter file
612 : ! in the paper this is Kab for B-Si
613 2008 : kab = 1.01_dp
614 : END SELECT
615 7308 : ELSEIF (za > 20 .AND. za < 30) THEN
616 : ! 3d
617 96 : IF (zb > 20 .AND. zb < 30) THEN
618 : ! 3d
619 : kab = 1.10_dp
620 76 : ELSEIF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
621 : ! 4d/5d/4f
622 0 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
623 : END IF
624 7212 : ELSEIF ((za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
625 : ! 4d/5d/4f
626 30 : IF (zb > 20 .AND. zb < 30) THEN
627 : ! 3d
628 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
629 30 : ELSEIF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
630 : ! 4d/5d/4f
631 28 : kab = 1.20_dp
632 : END IF
633 : END IF
634 : END IF
635 :
636 21816 : END FUNCTION xtb_set_kab
637 :
638 : END MODULE xtb_parameters
639 :
|
