(git:744416f)
Loading...
Searching...
No Matches
xtb_parameters.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Read xTB parameters.
10!> \author JGH (10.2018)
11! **************************************************************************************************
13
26 USE kinds, ONLY: default_string_length,&
27 dp
30 ptable
31 USE physcon, ONLY: bohr,&
32 evolt
34 USE xtb_types, ONLY: xtb_atom_type
35#include "./base/base_uses.f90"
36
37 IMPLICIT NONE
38
39 PRIVATE
40
41 INTEGER, PARAMETER, PRIVATE :: nelem = 106
42 ! H He
43 ! Li Be B C N O F Ne
44 ! Na Mg Al Si P S Cl Ar
45 ! K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
46 ! Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe
47 ! Cs Ba La Ce-Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
48 ! Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Ha 106
49
50!&<
51 ! Element Valence
52 INTEGER, DIMENSION(0:nelem), &
53 PARAMETER, PRIVATE :: zval = [-1, & ! 0
54 1, 2, & ! 2
55 1, 2, 3, 4, 5, 6, 7, 8, & ! 10
56 1, 2, 3, 4, 5, 6, 7, 8, & ! 18
57 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 36
58 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 54
59 1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, &
60 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 86
61 -1, -1, -1, 4, -1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]
62!&>
63
64!&<
65 ! Element Pauling Electronegativity
66 REAL(KIND=dp), DIMENSION(0:nelem), &
67 PARAMETER, PRIVATE :: eneg = [0.00_dp, & ! 0
68 2.20_dp, 3.00_dp, & ! 2
69 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
70 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
71 0.82_dp, 1.00_dp, 1.36_dp, 1.54_dp, 1.63_dp, 1.66_dp, 1.55_dp, 1.83_dp, &
72 1.88_dp, 1.91_dp, 1.90_dp, 1.65_dp, 1.81_dp, 2.01_dp, 2.18_dp, 2.55_dp, &
73 2.96_dp, 3.00_dp, & ! 36
74 0.82_dp, 0.95_dp, 1.22_dp, 1.33_dp, 1.60_dp, 2.16_dp, 1.90_dp, 2.20_dp, &
75 2.28_dp, 2.20_dp, 1.93_dp, 1.69_dp, 1.78_dp, 1.96_dp, 2.05_dp, 2.10_dp, &
76 2.66_dp, 2.60_dp, & ! 54
77 0.79_dp, 0.89_dp, 1.10_dp, &
78 1.12_dp, 1.13_dp, 1.14_dp, 1.15_dp, 1.17_dp, 1.18_dp, 1.20_dp, 1.21_dp, &
79 1.22_dp, 1.23_dp, 1.24_dp, 1.25_dp, 1.26_dp, 1.27_dp, & ! Lanthanides
80 1.30_dp, 1.50_dp, 2.36_dp, 1.90_dp, 2.20_dp, 2.20_dp, 2.28_dp, 2.54_dp, &
81 2.00_dp, 2.04_dp, 2.33_dp, 2.02_dp, 2.00_dp, 2.20_dp, 2.20_dp, & ! 86
82 0.70_dp, 0.89_dp, 1.10_dp, &
83 1.30_dp, 1.50_dp, 1.38_dp, 1.36_dp, 1.28_dp, 1.30_dp, 1.30_dp, 1.30_dp, &
84 1.30_dp, 1.30_dp, 1.30_dp, 1.30_dp, 1.30_dp, 1.50_dp, & ! Actinides
85 1.50_dp, 1.50_dp, 1.50_dp]
86!&>
87
88!&<
89 ! Shell occupation
90 INTEGER, DIMENSION(1:5, 0:nelem) :: occupation = reshape([0,0,0,0,0, & ! 0
91 1,0,0,0,0, 2,0,0,0,0, & ! 2
92 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
93 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
94 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, &
95 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
96 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, & !
97 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
98 1,0,0,0,0, 2,0,0,0,0, 2,0,1,0,0, &
99 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, &
100 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
101 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, &
102 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)
103 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, & !
104 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, &
105 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
106 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0], [5, nelem+1])
107!&>
108
109!&<
110 ! COVALENT RADII
111 ! based on "Atomic Radii of the Elements," M. Mantina, R. Valero, C. J. Cramer, and D. G. Truhlar,
112 ! in CRC Handbook of Chemistry and Physics, 91st Edition (2010-2011),
113 ! edited by W. M. Haynes (CRC Press, Boca Raton, FL, 2010), pages 9-49-9-50;
114 ! corrected Nov. 17, 2010 for the 92nd edition.
115 REAL(KIND=dp), DIMENSION(0:nelem), &
116 PARAMETER, PRIVATE :: crad = [0.00_dp, & ! 0
117 0.32_dp, 0.37_dp, & ! 2
118 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
119 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
120 2.00_dp, 1.74_dp, 1.59_dp, 1.48_dp, 1.44_dp, 1.30_dp, 1.29_dp, 1.24_dp, &
121 1.18_dp, 1.17_dp, 1.22_dp, 1.20_dp, 1.23_dp, 1.20_dp, 1.20_dp, 1.18_dp, &
122 1.17_dp, 1.16_dp, & ! 36
123 2.15_dp, 1.90_dp, 1.76_dp, 1.64_dp, 1.56_dp, 1.46_dp, 1.38_dp, 1.36_dp, &
124 1.34_dp, 1.30_dp, 1.36_dp, 1.40_dp, 1.42_dp, 1.40_dp, 1.40_dp, 1.37_dp, &
125 1.36_dp, 1.36_dp, & ! 54
126 2.38_dp, 2.06_dp, 1.94_dp, &
127 1.84_dp, 1.90_dp, 1.88_dp, 1.86_dp, 1.85_dp, 1.83_dp, 1.82_dp, 1.81_dp, &
128 1.80_dp, 1.79_dp, 1.77_dp, 1.77_dp, 1.78_dp, 1.74_dp, & ! Lanthanides
129 1.64_dp, 1.58_dp, 1.50_dp, 1.41_dp, 1.36_dp, 1.32_dp, 1.30_dp, 1.30_dp, &
130 1.32_dp, 1.44_dp, 1.45_dp, 1.50_dp, 1.42_dp, 1.48_dp, 1.46_dp, & ! 86
131 2.42_dp, 2.11_dp, 2.01_dp, &
132 1.90_dp, 1.84_dp, 1.83_dp, 1.80_dp, 1.80_dp, 1.51_dp, 0.96_dp, 1.54_dp, &
133 1.83_dp, 1.50_dp, 1.50_dp, 1.50_dp, 1.50_dp, 1.50_dp, & ! Actinides
134 1.50_dp, 1.50_dp, 1.50_dp]
135!&>
136
137!&<
138 ! Charge Limits (Mulliken)
139 REAL(KIND=dp), DIMENSION(0:nelem), &
140 PARAMETER, PRIVATE :: clmt = [0.00_dp, & ! 0
141 1.05_dp, 1.25_dp, & ! 2
142 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
143 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
144 1.05_dp, 2.05_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
145 3.50_dp, 3.50_dp, 3.50_dp, 2.50_dp, 2.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
146 1.25_dp, 1.00_dp, & ! 36
147 1.05_dp, 2.05_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
148 3.50_dp, 3.50_dp, 3.50_dp, 2.50_dp, 2.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
149 1.25_dp, 1.00_dp, & ! 54
150 1.05_dp, 2.05_dp, 3.00_dp, &
151 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, &
152 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, & ! Lanthanides
153 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
154 2.50_dp, 2.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 1.25_dp, 1.00_dp, & ! 86
155 1.05_dp, 2.05_dp, 3.00_dp, &
156 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, &
157 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, & ! Actinides
158 3.00_dp, 3.00_dp, 3.00_dp]
159!&>
160
161!&<
162 ! number of primitive gaussians per shell
163 INTEGER, PARAMETER :: number_of_primitives(1:3, 1:nelem) = reshape([&
164 & 4, 3, 0, 4, 0, 0, 6, 6, 0, 6, 6, 0, 6, 6, 0, 6, 6, 0, 6, 6, 0, &
165 & 6, 6, 0, 6, 6, 0, 6, 6, 4, 6, 6, 0, 6, 6, 0, 6, 6, 4, 6, 6, 4, &
166 & 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 0, 6, 6, 4, 4, 6, 6, &
167 & 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, &
168 & 4, 6, 6, 6, 6, 0, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, &
169 & 6, 6, 4, 6, 6, 0, 6, 6, 4, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, &
170 & 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 6, 6, 0, 6, 6, 4, &
171 & 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 0, 6, 6, 4, &
172 & 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, &
173 & 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, &
174 & 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, 4, 6, 6, &
175 & 4, 6, 6, 4, 6, 6, 6, 6, 0, 6, 6, 0, 6, 6, 0, 6, 6, 0, 6, 6, 4, &
176 & 6, 6, 4, 6, 6, 4, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, &
177 & 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, &
178 & 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, &
179 & 0, 0, 0], &
180 [3, nelem])
181!&>
182
183! *** Global parameters ***
184
185 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_parameters'
186
187! *** Public data types ***
188
191 PUBLIC :: metal, early3d, pp_gfn0
192
193CONTAINS
194
195! **************************************************************************************************
196!> \brief ...
197!> \param param ...
198!> \param gfn_type ...
199!> \param element_symbol ...
200!> \param parameter_file_path ...
201!> \param parameter_file_name ...
202!> \param para_env ...
203! **************************************************************************************************
204 SUBROUTINE xtb_parameters_init(param, gfn_type, element_symbol, &
205 parameter_file_path, parameter_file_name, &
206 para_env)
207
208 TYPE(xtb_atom_type), POINTER :: param
209 INTEGER, INTENT(IN) :: gfn_type
210 CHARACTER(LEN=2), INTENT(IN) :: element_symbol
211 CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
212 TYPE(mp_para_env_type), POINTER :: para_env
213
214 SELECT CASE (gfn_type)
215 CASE (0)
216 CALL xtb0_parameters_init(param, element_symbol, parameter_file_path, &
217 parameter_file_name, para_env)
218 CASE (1)
219 CALL xtb1_parameters_init(param, element_symbol, parameter_file_path, &
220 parameter_file_name, para_env)
221 CASE (2)
222 cpabort("gfn_type = 2 not yet supported")
223 CASE DEFAULT
224 cpabort("Wrong gfn_type")
225 END SELECT
226
227 END SUBROUTINE xtb_parameters_init
228
229! **************************************************************************************************
230!> \brief ...
231!> \param param ...
232!> \param element_symbol ...
233!> \param parameter_file_path ...
234!> \param parameter_file_name ...
235!> \param para_env ...
236! **************************************************************************************************
237 SUBROUTINE xtb0_parameters_init(param, element_symbol, parameter_file_path, parameter_file_name, &
238 para_env)
239
240 TYPE(xtb_atom_type), POINTER :: param
241 CHARACTER(LEN=2), INTENT(IN) :: element_symbol
242 CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
243 TYPE(mp_para_env_type), POINTER :: para_env
244
245 CHARACTER(len=2) :: esym
246 CHARACTER(len=default_string_length) :: aname, atag, filename
247 INTEGER :: i, l, zin, znum
248 LOGICAL :: at_end, found
249 TYPE(cp_parser_type) :: parser
250
251 filename = adjustl(trim(parameter_file_path))//adjustl(trim(parameter_file_name))
252 CALL parser_create(parser, filename, apply_preprocessing=.false., para_env=para_env)
253 found = .false.
254 znum = 0
255 CALL get_ptable_info(element_symbol, znum)
256 DO
257 at_end = .false.
258 CALL parser_get_next_line(parser, 1, at_end)
259 IF (at_end) EXIT
260 CALL parser_get_object(parser, aname)
261 CALL uppercase(aname)
262 IF (aname == "$Z") THEN
263 CALL parser_get_object(parser, zin)
264 IF (zin == znum) THEN
265 found = .true.
266 DO
267 CALL parser_get_next_line(parser, 1, at_end)
268 IF (at_end) THEN
269 cpabort("Incomplete xTB parameter file")
270 END IF
271 CALL parser_get_object(parser, aname)
272 CALL uppercase(aname)
273 SELECT CASE (aname)
274 CASE ("AO")
275 CALL parser_get_object(parser, atag)
276 CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
277 CASE ("LEV")
278 DO i = 1, param%nshell
279 CALL parser_get_object(parser, param%hen(i))
280 END DO
281 CASE ("EXP")
282 DO i = 1, param%nshell
283 CALL parser_get_object(parser, param%zeta(i))
284 END DO
285 CASE ("EN")
286 CALL parser_get_object(parser, param%en)
287 CASE ("GAM")
288 CALL parser_get_object(parser, param%eta)
289 CASE ("KQAT2")
290 CALL parser_get_object(parser, param%kqat2)
291 CASE ("KCNS")
292 CALL parser_get_object(parser, param%kcn(1))
293 param%kcn(1) = param%kcn(1)*0.1_dp !from orig xtb code
294 CASE ("KCNP")
295 CALL parser_get_object(parser, param%kcn(2))
296 param%kcn(2) = param%kcn(2)*0.1_dp !from orig xtb code
297 CASE ("KCND")
298 CALL parser_get_object(parser, param%kcn(3))
299 param%kcn(3) = param%kcn(3)*0.1_dp !from orig xtb code
300 CASE ("REPA")
301 CALL parser_get_object(parser, param%alpha)
302 CASE ("REPB")
303 CALL parser_get_object(parser, param%zneff)
304 CASE ("POLYS")
305 CALL parser_get_object(parser, param%kpoly(1))
306 CASE ("POLYP")
307 CALL parser_get_object(parser, param%kpoly(2))
308 CASE ("POLYD")
309 CALL parser_get_object(parser, param%kpoly(3))
310 CASE ("KQS")
311 CALL parser_get_object(parser, param%kq(1))
312 CASE ("KQP")
313 CALL parser_get_object(parser, param%kq(2))
314 CASE ("KQD")
315 CALL parser_get_object(parser, param%kq(3))
316 CASE ("XI")
317 CALL parser_get_object(parser, param%xi)
318 CASE ("KAPPA")
319 CALL parser_get_object(parser, param%kappa0)
320 CASE ("ALPG")
321 CALL parser_get_object(parser, param%alpg)
322 CASE ("$END")
323 EXIT
324 CASE DEFAULT
325 cpabort("Unknown parameter in xTB file")
326 END SELECT
327 END DO
328 ELSE
329 cycle
330 END IF
331 EXIT
332 END IF
333 END DO
334 IF (found) THEN
335 param%typ = "STANDARD"
336 param%symbol = element_symbol
337 param%defined = .true.
338 param%z = znum
339 param%aname = ptable(znum)%name
340 param%lmax = maxval(param%lval(1:param%nshell))
341 param%natorb = 0
342 DO i = 1, param%nshell
343 l = param%lval(i)
344 param%natorb = param%natorb + (2*l + 1)
345 END DO
346 param%zeff = zval(znum)
347 ELSE
348 esym = element_symbol
349 CALL uppercase(esym)
350 IF ("X " == esym) THEN
351 param%typ = "GHOST"
352 param%symbol = element_symbol
353 param%defined = .false.
354 param%z = 0
355 param%aname = "X "
356 param%lmax = 0
357 param%natorb = 0
358 param%nshell = 0
359 param%zeff = 0.0_dp
360 ELSE
361 param%defined = .false.
362 CALL cp_warn(__location__, "xTB parameters for element "//element_symbol// &
363 " were not found in the parameter file "//adjustl(trim(filename)))
364 END IF
365 END IF
366 CALL parser_release(parser)
367
368 END SUBROUTINE xtb0_parameters_init
369
370! **************************************************************************************************
371!> \brief ...
372!> \param param ...
373!> \param element_symbol ...
374!> \param parameter_file_path ...
375!> \param parameter_file_name ...
376!> \param para_env ...
377! **************************************************************************************************
378 SUBROUTINE xtb1_parameters_init(param, element_symbol, parameter_file_path, parameter_file_name, &
379 para_env)
380
381 TYPE(xtb_atom_type), POINTER :: param
382 CHARACTER(LEN=2), INTENT(IN) :: element_symbol
383 CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
384 TYPE(mp_para_env_type), POINTER :: para_env
385
386 CHARACTER(len=2) :: esym
387 CHARACTER(len=default_string_length) :: aname, atag, filename
388 INTEGER :: i, l, zin, znum
389 LOGICAL :: at_end, found
390 TYPE(cp_parser_type) :: parser
391
392 filename = adjustl(trim(parameter_file_path))//adjustl(trim(parameter_file_name))
393 CALL parser_create(parser, filename, apply_preprocessing=.false., para_env=para_env)
394 found = .false.
395 znum = 0
396 CALL get_ptable_info(element_symbol, znum)
397 DO
398 at_end = .false.
399 CALL parser_get_next_line(parser, 1, at_end)
400 IF (at_end) EXIT
401 CALL parser_get_object(parser, aname)
402 CALL uppercase(aname)
403 IF (aname == "$Z") THEN
404 CALL parser_get_object(parser, zin)
405 IF (zin == znum) THEN
406 found = .true.
407 DO
408 CALL parser_get_next_line(parser, 1, at_end)
409 IF (at_end) THEN
410 cpabort("Incomplete xTB parameter file")
411 END IF
412 CALL parser_get_object(parser, aname)
413 CALL uppercase(aname)
414 SELECT CASE (aname)
415 CASE ("AO")
416 CALL parser_get_object(parser, atag)
417 CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
418 CASE ("LEV")
419 DO i = 1, param%nshell
420 CALL parser_get_object(parser, param%hen(i))
421 END DO
422 CASE ("EXP")
423 DO i = 1, param%nshell
424 CALL parser_get_object(parser, param%zeta(i))
425 END DO
426 CASE ("GAM")
427 CALL parser_get_object(parser, param%eta)
428 CASE ("GAM3")
429 CALL parser_get_object(parser, param%xgamma)
430 CASE ("CXB")
431 CALL parser_get_object(parser, param%kx)
432 CASE ("REPA")
433 CALL parser_get_object(parser, param%alpha)
434 CASE ("REPB")
435 CALL parser_get_object(parser, param%zneff)
436 CASE ("POLYS")
437 CALL parser_get_object(parser, param%kpoly(1))
438 CASE ("POLYP")
439 CALL parser_get_object(parser, param%kpoly(2))
440 CASE ("POLYD")
441 CALL parser_get_object(parser, param%kpoly(3))
442 CASE ("LPARP")
443 CALL parser_get_object(parser, param%kappa(2))
444 CASE ("LPARD")
445 CALL parser_get_object(parser, param%kappa(3))
446 CASE ("$END")
447 EXIT
448 CASE DEFAULT
449 cpabort("Unknown parameter in xTB file")
450 END SELECT
451 END DO
452 ELSE
453 cycle
454 END IF
455 EXIT
456 END IF
457 END DO
458 IF (found) THEN
459 param%typ = "STANDARD"
460 param%symbol = element_symbol
461 param%defined = .true.
462 param%z = znum
463 param%aname = ptable(znum)%name
464 param%lmax = maxval(param%lval(1:param%nshell))
465 param%natorb = 0
466 DO i = 1, param%nshell
467 l = param%lval(i)
468 param%natorb = param%natorb + (2*l + 1)
469 END DO
470 param%zeff = zval(znum)
471 ELSE
472 esym = element_symbol
473 CALL uppercase(esym)
474 IF ("X " == esym) THEN
475 param%typ = "GHOST"
476 param%symbol = element_symbol
477 param%defined = .false.
478 param%z = 0
479 param%aname = "X "
480 param%lmax = 0
481 param%natorb = 0
482 param%nshell = 0
483 param%zeff = 0.0_dp
484 ELSE
485 param%defined = .false.
486 CALL cp_warn(__location__, "xTB parameters for element "//element_symbol// &
487 " were not found in the parameter file "//adjustl(trim(filename)))
488 END IF
489 END IF
490 CALL parser_release(parser)
491
492 END SUBROUTINE xtb1_parameters_init
493
494! **************************************************************************************************
495!> \brief ...
496!> \param param ...
497!> \param gfn_type ...
498!> \param element_symbol ...
499!> \param parameter_file_path ...
500!> \param spinpol_param_file_name ...
501!> \param para_env ...
502! **************************************************************************************************
503 SUBROUTINE xtb_spinpol_init(param, gfn_type, element_symbol, parameter_file_path, spinpol_param_file_name, &
504 para_env)
505
506 TYPE(xtb_atom_type), POINTER :: param
507 INTEGER, INTENT(IN) :: gfn_type
508 CHARACTER(LEN=2), INTENT(IN) :: element_symbol
509 CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, &
510 spinpol_param_file_name
511 TYPE(mp_para_env_type), POINTER :: para_env
512
513 CHARACTER(len=default_string_length) :: filename
514 INTEGER :: zin, znum
515 LOGICAL :: at_end
516 TYPE(cp_parser_type) :: parser
517
518 SELECT CASE (gfn_type)
519 CASE (0)
520 cpabort("gfn_type = 0: No spin polarisation possible!")
521 CASE (1)
522 ! OK
523 CASE (2)
524 cpabort("gfn_type = 2 not yet supported")
525 CASE DEFAULT
526 cpabort("Wrong gfn_type")
527 END SELECT
528
529 filename = adjustl(trim(parameter_file_path))//adjustl(trim(spinpol_param_file_name))
530 CALL parser_create(parser, filename, apply_preprocessing=.false., para_env=para_env)
531 znum = 0
532 param%wall = 0.0_dp
533 CALL get_ptable_info(element_symbol, znum)
534 DO
535 at_end = .false.
536 CALL parser_get_next_line(parser, 1, at_end)
537 IF (at_end) EXIT
538 CALL parser_get_object(parser, zin)
539 IF (zin == znum) THEN
540 CALL parser_get_object(parser, param%wall(1, 1))
541 CALL parser_get_object(parser, param%wall(1, 2))
542 CALL parser_get_object(parser, param%wall(2, 2))
543 CALL parser_get_object(parser, param%wall(1, 3))
544 CALL parser_get_object(parser, param%wall(2, 3))
545 CALL parser_get_object(parser, param%wall(3, 3))
546 param%wall(2, 1) = param%wall(1, 2)
547 param%wall(3, 1) = param%wall(1, 3)
548 param%wall(3, 2) = param%wall(2, 3)
549 END IF
550 END DO
551 CALL parser_release(parser)
552
553 END SUBROUTINE xtb_spinpol_init
554
555! **************************************************************************************************
556!> \brief ...
557!> \param param ...
558!> \param gfn_type ...
559!> \param xtb_control ...
560! **************************************************************************************************
561 SUBROUTINE xtb_spinpol_ext(param, gfn_type, xtb_control)
562 TYPE(xtb_atom_type), POINTER :: param
563 INTEGER, INTENT(IN) :: gfn_type
564 TYPE(xtb_control_type), INTENT(IN), POINTER :: xtb_control
565
566 INTEGER :: i
567
568 SELECT CASE (gfn_type)
569 CASE (0)
570 cpabort("gfn_type = 0: No spin polarisation possible!")
571 CASE (1)
572 ! OK
573 CASE (2)
574 cpabort("gfn_type = 2 not yet supported")
575 CASE DEFAULT
576 cpabort("Wrong gfn_type")
577 END SELECT
578
579 IF (param%defined) THEN
580 IF (ASSOCIATED(xtb_control%spinpol_type)) THEN
581 DO i = 1, SIZE(xtb_control%spinpol_type)
582 IF (xtb_control%spinpol_type(i) == param%z) THEN
583 param%wall(1, 1) = xtb_control%spinpol_vals(1, i)
584 param%wall(1, 2) = xtb_control%spinpol_vals(2, i)
585 param%wall(2, 2) = xtb_control%spinpol_vals(3, i)
586 param%wall(1, 3) = xtb_control%spinpol_vals(4, i)
587 param%wall(2, 3) = xtb_control%spinpol_vals(5, i)
588 param%wall(3, 3) = xtb_control%spinpol_vals(6, i)
589 param%wall(2, 1) = param%wall(1, 2)
590 param%wall(3, 1) = param%wall(1, 3)
591 param%wall(3, 2) = param%wall(2, 3)
592 EXIT
593 END IF
594 END DO
595 END IF
596 END IF
597
598 END SUBROUTINE xtb_spinpol_ext
599
600! **************************************************************************************************
601!> \brief Read atom parameters for xTB Hamiltonian from input file
602!> \param param ...
603! **************************************************************************************************
604 SUBROUTINE xtb_parameters_set(param)
605
606 TYPE(xtb_atom_type), POINTER :: param
607
608 INTEGER :: i, is, l, na
609 REAL(kind=dp), DIMENSION(5) :: kp
610
611 IF (param%defined) THEN
612 ! AO to shell pointer
613 ! AO to l-qn pointer
614 na = 0
615 DO is = 1, param%nshell
616 l = param%lval(is)
617 DO i = 1, 2*l + 1
618 na = na + 1
619 param%nao(na) = is
620 param%lao(na) = l
621 END DO
622 END DO
623 !
624 i = param%z
625 ! Electronegativity
626 param%electronegativity = eneg(i)
627 IF (param%en == 0.0_dp) param%en = eneg(i)
628 ! covalent radius
629 param%rcov = crad(i)*bohr
630 ! shell occupances
631 param%occupation(:) = occupation(:, i)
632 ! number of primitive Gaussians per shell
633 param%ngauss = 0
634 param%ngauss(1:3) = number_of_primitives(:, i)
635 ! check for consistency
636 IF (abs(param%zeff - sum(param%occupation)) > 1.e-10_dp) THEN
637 CALL cp_abort(__location__, "Element <"//trim(param%aname)//"> has inconsistent shell occupations")
638 END IF
639 ! orbital energies [evolt] -> [a.u.]
640 param%hen = param%hen/evolt
641 ! some forgotten scaling parameters (not in orig. paper)
642 param%xgamma = 0.1_dp*param%xgamma
643 param%kpoly(:) = 0.01_dp*param%kpoly(:)
644 param%kappa(:) = 0.1_dp*param%kappa(:)
645 ! we have 1/6 g * q**3 (not 1/3)
646 param%xgamma = -2.0_dp*param%xgamma
647 ! we need kpoly in shell order
648 kp(:) = param%kpoly(:)
649 param%kpoly(:) = 0.0_dp
650 DO is = 1, param%nshell
651 l = param%lval(is)
652 param%kpoly(is) = kp(l + 1)
653 END DO
654 ! kx
655 param%kx = 0.1_dp*param%kx
656 IF (param%kx < -5._dp) THEN
657 ! use defaults
658 SELECT CASE (param%z)
659 CASE DEFAULT
660 param%kx = 0.0_dp
661 CASE (35) ! Br
662 param%kx = 0.1_dp*0.381742_dp
663 CASE (53) ! I
664 param%kx = 0.1_dp*0.321944_dp
665 CASE (85) ! At
666 param%kx = 0.1_dp*0.220000_dp
667 END SELECT
668 END IF
669 ! chmax
670 param%chmax = clmt(i)
671 END IF
672
673 END SUBROUTINE xtb_parameters_set
674
675! **************************************************************************************************
676!> \brief ...
677!> \param param ...
678!> \param gto_basis_set ...
679!> \param ngauss ...
680!> \param ngaussflex ...
681! **************************************************************************************************
682 SUBROUTINE init_xtb_basis(param, gto_basis_set, ngauss, ngaussflex)
683
684 TYPE(xtb_atom_type), POINTER :: param
685 TYPE(gto_basis_set_type), POINTER :: gto_basis_set
686 INTEGER, INTENT(IN) :: ngauss
687 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: ngaussflex
688
689 CHARACTER(LEN=6), DIMENSION(:), POINTER :: symbol
690 INTEGER :: i, nshell
691 INTEGER, DIMENSION(:), POINTER :: lq, nq
692 REAL(kind=dp), DIMENSION(:), POINTER :: zet
693 TYPE(sto_basis_set_type), POINTER :: sto_basis_set
694
695 IF (ASSOCIATED(param)) THEN
696 IF (param%defined) THEN
697 NULLIFY (sto_basis_set)
698 CALL allocate_sto_basis_set(sto_basis_set)
699 nshell = param%nshell
700
701 ALLOCATE (symbol(1:nshell))
702 symbol = ""
703 DO i = 1, nshell
704 SELECT CASE (param%lval(i))
705 CASE (0)
706 WRITE (symbol(i), '(I1,A1)') param%nval(i), "S"
707 CASE (1)
708 WRITE (symbol(i), '(I1,A1)') param%nval(i), "P"
709 CASE (2)
710 WRITE (symbol(i), '(I1,A1)') param%nval(i), "D"
711 CASE (3)
712 WRITE (symbol(i), '(I1,A1)') param%nval(i), "F"
713 CASE DEFAULT
714 cpabort('BASIS SET OUT OF RANGE (lval)')
715 END SELECT
716 END DO
717
718 IF (nshell > 0) THEN
719 ALLOCATE (nq(nshell), lq(nshell), zet(nshell))
720 nq(1:nshell) = param%nval(1:nshell)
721 lq(1:nshell) = param%lval(1:nshell)
722 zet(1:nshell) = param%zeta(1:nshell)
723 CALL set_sto_basis_set(sto_basis_set, name=param%aname, nshell=nshell, symbol=symbol, &
724 nq=nq, lq=lq, zet=zet)
725 IF (PRESENT(ngaussflex)) THEN
726 CALL create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss=ngauss, ortho=.true., &
727 ngaussflex=ngaussflex)
728 ELSE
729 CALL create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss=ngauss, ortho=.true.)
730 END IF
731 END IF
732
733 ! this will remove the allocated arrays
734 CALL deallocate_sto_basis_set(sto_basis_set)
735 DEALLOCATE (symbol, nq, lq, zet)
736 END IF
737
738 ELSE
739 cpabort("The pointer param is not associated")
740 END IF
741
742 END SUBROUTINE init_xtb_basis
743
744! **************************************************************************************************
745!> \brief ...
746!> \param za ...
747!> \param zb ...
748!> \param xtb_control ...
749!> \return ...
750! **************************************************************************************************
751 FUNCTION xtb_set_kab(za, zb, xtb_control) RESULT(kab)
752
753 INTEGER, INTENT(IN) :: za, zb
754 TYPE(xtb_control_type), INTENT(IN), POINTER :: xtb_control
755 REAL(kind=dp) :: kab
756
757 INTEGER :: j, z
758 LOGICAL :: custom
759
760 kab = 1.0_dp
761 custom = .false.
762
763 IF (xtb_control%kab_nval > 0) THEN
764 DO j = 1, xtb_control%kab_nval
765 IF ((za == xtb_control%kab_types(1, j) .AND. &
766 zb == xtb_control%kab_types(2, j)) .OR. &
767 (za == xtb_control%kab_types(2, j) .AND. &
768 zb == xtb_control%kab_types(1, j))) THEN
769 custom = .true.
770 kab = xtb_control%kab_vals(j)
771 EXIT
772 END IF
773 END DO
774 END IF
775
776 IF (.NOT. custom) THEN
777 IF (za == 1 .OR. zb == 1) THEN
778 ! hydrogen
779 z = za + zb - 1
780 SELECT CASE (z)
781 CASE (1)
782 kab = 0.96_dp
783 CASE (5)
784 kab = 0.95_dp
785 CASE (7)
786 kab = 1.04_dp
787 CASE (28)
788 kab = 0.90_dp
789 CASE (75)
790 kab = 0.80_dp
791 CASE (78)
792 kab = 0.80_dp
793 END SELECT
794 ELSE IF (za == 5 .OR. zb == 5) THEN
795 ! Boron
796 z = za + zb - 5
797 SELECT CASE (z)
798 CASE (15)
799 kab = 0.97_dp
800 END SELECT
801 ELSE IF (za == 7 .OR. zb == 7) THEN
802 ! Nitrogen
803 z = za + zb - 7
804 SELECT CASE (z)
805 CASE (14)
806 !xtb orig code parameter file
807 ! in the paper this is Kab for B-Si
808 kab = 1.01_dp
809 END SELECT
810 ELSE IF (za > 20 .AND. za < 30) THEN
811 ! 3d
812 IF (zb > 20 .AND. zb < 30) THEN
813 ! 3d
814 kab = 1.10_dp
815 ELSE IF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
816 ! 4d/5d/4f
817 kab = 0.50_dp*(1.20_dp + 1.10_dp)
818 END IF
819 ELSE IF ((za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
820 ! 4d/5d/4f
821 IF (zb > 20 .AND. zb < 30) THEN
822 ! 3d
823 kab = 0.50_dp*(1.20_dp + 1.10_dp)
824 ELSE IF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
825 ! 4d/5d/4f
826 kab = 1.20_dp
827 END IF
828 END IF
829 END IF
830
831 END FUNCTION xtb_set_kab
832
833! **************************************************************************************************
834!> \brief ...
835!> \param atag ...
836!> \param nshell ...
837!> \param nval ...
838!> \param lval ...
839!> \return ...
840! **************************************************************************************************
841 SUBROUTINE xtb_get_shells(atag, nshell, nval, lval)
842 CHARACTER(len=*) :: atag
843 INTEGER :: nshell
844 INTEGER, DIMENSION(:) :: nval, lval
845
846 CHARACTER(LEN=1) :: ltag
847 CHARACTER(LEN=10) :: aotag
848 INTEGER :: i, j
849
850 aotag = adjustl(trim(atag))
851 nshell = len(trim(aotag))/2
852 DO i = 1, nshell
853 j = (i - 1)*2 + 1
854 READ (aotag(j:j), fmt="(i1)") nval(i)
855 READ (aotag(j + 1:j + 1), fmt="(A1)") ltag
856 CALL uppercase(ltag)
857 SELECT CASE (ltag)
858 CASE ("S")
859 lval(i) = 0
860 CASE ("P")
861 lval(i) = 1
862 CASE ("D")
863 lval(i) = 2
864 CASE DEFAULT
865 END SELECT
866 END DO
867
868 END SUBROUTINE xtb_get_shells
869
870! **************************************************************************************************
871!> \brief ...
872!> \param z ...
873!> \return ...
874! **************************************************************************************************
875 FUNCTION metal(z) RESULT(ismetal)
876 INTEGER :: z
877 LOGICAL :: ismetal
878
879 SELECT CASE (z)
880 CASE DEFAULT
881 ismetal = .true.
882 CASE (1:2, 6:10, 14:18, 32:36, 50:54, 82:86)
883 ismetal = .false.
884 END SELECT
885
886 END FUNCTION metal
887
888! **************************************************************************************************
889!> \brief ...
890!> \param z ...
891!> \return ...
892! **************************************************************************************************
893 FUNCTION early3d(z) RESULT(isearly3d)
894 INTEGER :: z
895 LOGICAL :: isearly3d
896
897 isearly3d = .false.
898 IF (z >= 21 .AND. z <= 24) isearly3d = .true.
899
900 END FUNCTION early3d
901
902! **************************************************************************************************
903!> \brief ...
904!> \param za ...
905!> \param zb ...
906!> \return ...
907! **************************************************************************************************
908 FUNCTION pp_gfn0(za, zb) RESULT(pparm)
909 INTEGER :: za, zb
910 REAL(kind=dp) :: pparm
911
912 pparm = 1.0_dp
913 IF ((za > 20 .AND. za < 30) .OR. (za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
914 IF ((zb > 20 .AND. zb < 30) .OR. (zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
915 pparm = 1.1_dp
916 IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
917 IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
918 pparm = 0.9_dp
919 END IF
920 END IF
921 END IF
922 END IF
923
924 END FUNCTION pp_gfn0
925
926END MODULE xtb_parameters
subroutine, public create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho, ngaussflex)
...
subroutine, public deallocate_sto_basis_set(sto_basis_set)
...
subroutine, public allocate_sto_basis_set(sto_basis_set)
...
subroutine, public set_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet)
...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Interface to the message passing library MPI.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
real(kind=dp), parameter, public bohr
Definition physcon.F:147
Utilities for string manipulations.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Read xTB parameters.
subroutine, public xtb_parameters_set(param)
Read atom parameters for xTB Hamiltonian from input file.
logical function, public metal(z)
...
logical function, public early3d(z)
...
subroutine, public xtb_parameters_init(param, gfn_type, element_symbol, parameter_file_path, parameter_file_name, para_env)
...
subroutine, public xtb_spinpol_ext(param, gfn_type, xtb_control)
...
real(kind=dp) function, public pp_gfn0(za, zb)
...
subroutine, public xtb_spinpol_init(param, gfn_type, element_symbol, parameter_file_path, spinpol_param_file_name, para_env)
...
subroutine, public init_xtb_basis(param, gto_basis_set, ngauss, ngaussflex)
...
real(kind=dp) function, public xtb_set_kab(za, zb, xtb_control)
...
Definition of the xTB parameter types.
Definition xtb_types.F:20
stores all the informations relevant to an mpi environment