Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_parser_methods, ONLY: parser_get_next_line,&
22 : parser_get_object
23 : USE cp_parser_types, ONLY: cp_parser_type,&
24 : parser_create,&
25 : parser_release
26 : USE kinds, ONLY: default_string_length,&
27 : dp
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE periodic_table, ONLY: get_ptable_info,&
30 : ptable
31 : USE physcon, ONLY: bohr,&
32 : evolt
33 : USE string_utilities, ONLY: uppercase
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 : ! *** Global parameters ***
162 :
163 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_parameters'
164 :
165 : ! *** Public data types ***
166 :
167 : PUBLIC :: xtb_parameters_init, xtb_parameters_set, init_xtb_basis, xtb_set_kab
168 : PUBLIC :: xtb_spinpol_init, xtb_spinpol_ext
169 : PUBLIC :: metal, early3d, pp_gfn0
170 :
171 : CONTAINS
172 :
173 : ! **************************************************************************************************
174 : !> \brief ...
175 : !> \param param ...
176 : !> \param gfn_type ...
177 : !> \param element_symbol ...
178 : !> \param parameter_file_path ...
179 : !> \param parameter_file_name ...
180 : !> \param para_env ...
181 : ! **************************************************************************************************
182 2256 : SUBROUTINE xtb_parameters_init(param, gfn_type, element_symbol, &
183 : parameter_file_path, parameter_file_name, &
184 : para_env)
185 :
186 : TYPE(xtb_atom_type), POINTER :: param
187 : INTEGER, INTENT(IN) :: gfn_type
188 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
189 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
190 : TYPE(mp_para_env_type), POINTER :: para_env
191 :
192 3730 : SELECT CASE (gfn_type)
193 : CASE (0)
194 : CALL xtb0_parameters_init(param, element_symbol, parameter_file_path, &
195 1474 : parameter_file_name, para_env)
196 : CASE (1)
197 : CALL xtb1_parameters_init(param, element_symbol, parameter_file_path, &
198 782 : parameter_file_name, para_env)
199 : CASE (2)
200 0 : CPABORT("gfn_type = 2 not yet supported")
201 : CASE DEFAULT
202 2256 : CPABORT("Wrong gfn_type")
203 : END SELECT
204 :
205 2256 : END SUBROUTINE xtb_parameters_init
206 :
207 : ! **************************************************************************************************
208 : !> \brief ...
209 : !> \param param ...
210 : !> \param element_symbol ...
211 : !> \param parameter_file_path ...
212 : !> \param parameter_file_name ...
213 : !> \param para_env ...
214 : ! **************************************************************************************************
215 2948 : SUBROUTINE xtb0_parameters_init(param, element_symbol, parameter_file_path, parameter_file_name, &
216 : para_env)
217 :
218 : TYPE(xtb_atom_type), POINTER :: param
219 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
220 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
221 : TYPE(mp_para_env_type), POINTER :: para_env
222 :
223 : CHARACTER(len=2) :: esym
224 : CHARACTER(len=default_string_length) :: aname, atag, filename
225 : INTEGER :: i, l, zin, znum
226 : LOGICAL :: at_end, found
227 : TYPE(cp_parser_type) :: parser
228 :
229 1474 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(parameter_file_name))
230 1474 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
231 1474 : found = .FALSE.
232 : znum = 0
233 1474 : CALL get_ptable_info(element_symbol, znum)
234 : DO
235 : at_end = .FALSE.
236 554742 : CALL parser_get_next_line(parser, 1, at_end)
237 554742 : IF (at_end) EXIT
238 554742 : CALL parser_get_object(parser, aname)
239 554742 : CALL uppercase(aname)
240 554742 : IF (aname == "$Z") THEN
241 27298 : CALL parser_get_object(parser, zin)
242 27298 : IF (zin == znum) THEN
243 26628 : found = .TRUE.
244 : DO
245 26628 : CALL parser_get_next_line(parser, 1, at_end)
246 26628 : IF (at_end) THEN
247 0 : CPABORT("Incomplete xTB parameter file")
248 : END IF
249 26628 : CALL parser_get_object(parser, aname)
250 26628 : CALL uppercase(aname)
251 1474 : SELECT CASE (aname)
252 : CASE ("AO")
253 1474 : CALL parser_get_object(parser, atag)
254 1474 : CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
255 : CASE ("LEV")
256 4960 : DO i = 1, param%nshell
257 4960 : CALL parser_get_object(parser, param%hen(i))
258 : END DO
259 : CASE ("EXP")
260 4960 : DO i = 1, param%nshell
261 4960 : CALL parser_get_object(parser, param%zeta(i))
262 : END DO
263 : CASE ("EN")
264 616 : CALL parser_get_object(parser, param%en)
265 : CASE ("GAM")
266 1474 : CALL parser_get_object(parser, param%eta)
267 : CASE ("KQAT2")
268 1474 : CALL parser_get_object(parser, param%kqat2)
269 : CASE ("KCNS")
270 1474 : CALL parser_get_object(parser, param%kcn(1))
271 1474 : param%kcn(1) = param%kcn(1)*0.1_dp !from orig xtb code
272 : CASE ("KCNP")
273 1254 : CALL parser_get_object(parser, param%kcn(2))
274 1254 : param%kcn(2) = param%kcn(2)*0.1_dp !from orig xtb code
275 : CASE ("KCND")
276 538 : CALL parser_get_object(parser, param%kcn(3))
277 538 : param%kcn(3) = param%kcn(3)*0.1_dp !from orig xtb code
278 : CASE ("REPA")
279 1474 : CALL parser_get_object(parser, param%alpha)
280 : CASE ("REPB")
281 1474 : CALL parser_get_object(parser, param%zneff)
282 : CASE ("POLYS")
283 1474 : CALL parser_get_object(parser, param%kpoly(1))
284 : CASE ("POLYP")
285 1254 : CALL parser_get_object(parser, param%kpoly(2))
286 : CASE ("POLYD")
287 538 : CALL parser_get_object(parser, param%kpoly(3))
288 : CASE ("KQS")
289 1474 : CALL parser_get_object(parser, param%kq(1))
290 : CASE ("KQP")
291 1254 : CALL parser_get_object(parser, param%kq(2))
292 : CASE ("KQD")
293 538 : CALL parser_get_object(parser, param%kq(3))
294 : CASE ("XI")
295 1474 : CALL parser_get_object(parser, param%xi)
296 : CASE ("KAPPA")
297 1474 : CALL parser_get_object(parser, param%kappa0)
298 : CASE ("ALPG")
299 1474 : CALL parser_get_object(parser, param%alpg)
300 : CASE ("$END")
301 0 : EXIT
302 : CASE DEFAULT
303 26628 : CPABORT("Unknown parameter in xTB file")
304 : END SELECT
305 : END DO
306 : ELSE
307 : CYCLE
308 : END IF
309 : EXIT
310 : END IF
311 : END DO
312 1474 : IF (found) THEN
313 1474 : param%typ = "STANDARD"
314 1474 : param%symbol = element_symbol
315 1474 : param%defined = .TRUE.
316 1474 : param%z = znum
317 1474 : param%aname = ptable(znum)%name
318 4960 : param%lmax = MAXVAL(param%lval(1:param%nshell))
319 1474 : param%natorb = 0
320 4960 : DO i = 1, param%nshell
321 3486 : l = param%lval(i)
322 4960 : param%natorb = param%natorb + (2*l + 1)
323 : END DO
324 1474 : param%zeff = zval(znum)
325 : ELSE
326 0 : esym = element_symbol
327 0 : CALL uppercase(esym)
328 0 : IF ("X " == esym) THEN
329 0 : param%typ = "GHOST"
330 0 : param%symbol = element_symbol
331 0 : param%defined = .FALSE.
332 0 : param%z = 0
333 0 : param%aname = "X "
334 0 : param%lmax = 0
335 0 : param%natorb = 0
336 0 : param%nshell = 0
337 0 : param%zeff = 0.0_dp
338 : ELSE
339 0 : param%defined = .FALSE.
340 : CALL cp_warn(__LOCATION__, "xTB parameters for element "//element_symbol// &
341 0 : " were not found in the parameter file "//ADJUSTL(TRIM(filename)))
342 : END IF
343 : END IF
344 1474 : CALL parser_release(parser)
345 :
346 4422 : END SUBROUTINE xtb0_parameters_init
347 :
348 : ! **************************************************************************************************
349 : !> \brief ...
350 : !> \param param ...
351 : !> \param element_symbol ...
352 : !> \param parameter_file_path ...
353 : !> \param parameter_file_name ...
354 : !> \param para_env ...
355 : ! **************************************************************************************************
356 1564 : SUBROUTINE xtb1_parameters_init(param, element_symbol, parameter_file_path, parameter_file_name, &
357 : para_env)
358 :
359 : TYPE(xtb_atom_type), POINTER :: param
360 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
361 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
362 : TYPE(mp_para_env_type), POINTER :: para_env
363 :
364 : CHARACTER(len=2) :: esym
365 : CHARACTER(len=default_string_length) :: aname, atag, filename
366 : INTEGER :: i, l, zin, znum
367 : LOGICAL :: at_end, found
368 : TYPE(cp_parser_type) :: parser
369 :
370 782 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(parameter_file_name))
371 782 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
372 782 : found = .FALSE.
373 : znum = 0
374 782 : CALL get_ptable_info(element_symbol, znum)
375 : DO
376 : at_end = .FALSE.
377 90874 : CALL parser_get_next_line(parser, 1, at_end)
378 90874 : IF (at_end) EXIT
379 90872 : CALL parser_get_object(parser, aname)
380 90872 : CALL uppercase(aname)
381 90874 : IF (aname == "$Z") THEN
382 6002 : CALL parser_get_object(parser, zin)
383 6002 : IF (zin == znum) THEN
384 7570 : found = .TRUE.
385 : DO
386 7570 : CALL parser_get_next_line(parser, 1, at_end)
387 7570 : IF (at_end) THEN
388 0 : CPABORT("Incomplete xTB parameter file")
389 : END IF
390 7570 : CALL parser_get_object(parser, aname)
391 7570 : CALL uppercase(aname)
392 780 : SELECT CASE (aname)
393 : CASE ("AO")
394 780 : CALL parser_get_object(parser, atag)
395 780 : CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
396 : CASE ("LEV")
397 2436 : DO i = 1, param%nshell
398 2436 : CALL parser_get_object(parser, param%hen(i))
399 : END DO
400 : CASE ("EXP")
401 2436 : DO i = 1, param%nshell
402 2436 : CALL parser_get_object(parser, param%zeta(i))
403 : END DO
404 : CASE ("GAM")
405 780 : CALL parser_get_object(parser, param%eta)
406 : CASE ("GAM3")
407 488 : CALL parser_get_object(parser, param%xgamma)
408 : CASE ("CXB")
409 12 : CALL parser_get_object(parser, param%kx)
410 : CASE ("REPA")
411 780 : CALL parser_get_object(parser, param%alpha)
412 : CASE ("REPB")
413 780 : CALL parser_get_object(parser, param%zneff)
414 : CASE ("POLYS")
415 488 : CALL parser_get_object(parser, param%kpoly(1))
416 : CASE ("POLYP")
417 488 : CALL parser_get_object(parser, param%kpoly(2))
418 : CASE ("POLYD")
419 96 : CALL parser_get_object(parser, param%kpoly(3))
420 : CASE ("LPARP")
421 478 : CALL parser_get_object(parser, param%kappa(2))
422 : CASE ("LPARD")
423 60 : CALL parser_get_object(parser, param%kappa(3))
424 : CASE ("$END")
425 0 : EXIT
426 : CASE DEFAULT
427 7570 : CPABORT("Unknown parameter in xTB file")
428 : END SELECT
429 : END DO
430 : ELSE
431 : CYCLE
432 : END IF
433 : EXIT
434 : END IF
435 : END DO
436 782 : IF (found) THEN
437 780 : param%typ = "STANDARD"
438 780 : param%symbol = element_symbol
439 780 : param%defined = .TRUE.
440 780 : param%z = znum
441 780 : param%aname = ptable(znum)%name
442 2436 : param%lmax = MAXVAL(param%lval(1:param%nshell))
443 780 : param%natorb = 0
444 2436 : DO i = 1, param%nshell
445 1656 : l = param%lval(i)
446 2436 : param%natorb = param%natorb + (2*l + 1)
447 : END DO
448 780 : param%zeff = zval(znum)
449 : ELSE
450 2 : esym = element_symbol
451 2 : CALL uppercase(esym)
452 2 : IF ("X " == esym) THEN
453 2 : param%typ = "GHOST"
454 2 : param%symbol = element_symbol
455 2 : param%defined = .FALSE.
456 2 : param%z = 0
457 2 : param%aname = "X "
458 2 : param%lmax = 0
459 2 : param%natorb = 0
460 2 : param%nshell = 0
461 2 : param%zeff = 0.0_dp
462 : ELSE
463 0 : param%defined = .FALSE.
464 : CALL cp_warn(__LOCATION__, "xTB parameters for element "//element_symbol// &
465 0 : " were not found in the parameter file "//ADJUSTL(TRIM(filename)))
466 : END IF
467 : END IF
468 782 : CALL parser_release(parser)
469 :
470 2346 : END SUBROUTINE xtb1_parameters_init
471 :
472 : ! **************************************************************************************************
473 : !> \brief ...
474 : !> \param param ...
475 : !> \param gfn_type ...
476 : !> \param element_symbol ...
477 : !> \param parameter_file_path ...
478 : !> \param spinpol_param_file_name ...
479 : !> \param para_env ...
480 : ! **************************************************************************************************
481 116 : SUBROUTINE xtb_spinpol_init(param, gfn_type, element_symbol, parameter_file_path, spinpol_param_file_name, &
482 : para_env)
483 :
484 : TYPE(xtb_atom_type), POINTER :: param
485 : INTEGER, INTENT(IN) :: gfn_type
486 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
487 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, &
488 : spinpol_param_file_name
489 : TYPE(mp_para_env_type), POINTER :: para_env
490 :
491 : CHARACTER(len=default_string_length) :: filename
492 : INTEGER :: zin, znum
493 : LOGICAL :: at_end
494 : TYPE(cp_parser_type) :: parser
495 :
496 58 : SELECT CASE (gfn_type)
497 : CASE (0)
498 0 : CPABORT("gfn_type = 0: No spin polarisation possible!")
499 : CASE (1)
500 : ! OK
501 : CASE (2)
502 0 : CPABORT("gfn_type = 2 not yet supported")
503 : CASE DEFAULT
504 0 : CPABORT("Wrong gfn_type")
505 : END SELECT
506 :
507 58 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(spinpol_param_file_name))
508 58 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
509 : znum = 0
510 754 : param%wall = 0.0_dp
511 58 : CALL get_ptable_info(element_symbol, znum)
512 4904 : DO
513 : at_end = .FALSE.
514 4962 : CALL parser_get_next_line(parser, 1, at_end)
515 4962 : IF (at_end) EXIT
516 4904 : CALL parser_get_object(parser, zin)
517 4962 : IF (zin == znum) THEN
518 58 : CALL parser_get_object(parser, param%wall(1, 1))
519 58 : CALL parser_get_object(parser, param%wall(1, 2))
520 58 : CALL parser_get_object(parser, param%wall(2, 2))
521 58 : CALL parser_get_object(parser, param%wall(1, 3))
522 58 : CALL parser_get_object(parser, param%wall(2, 3))
523 58 : CALL parser_get_object(parser, param%wall(3, 3))
524 58 : param%wall(2, 1) = param%wall(1, 2)
525 58 : param%wall(3, 1) = param%wall(1, 3)
526 58 : param%wall(3, 2) = param%wall(2, 3)
527 : END IF
528 : END DO
529 58 : CALL parser_release(parser)
530 :
531 174 : END SUBROUTINE xtb_spinpol_init
532 :
533 : ! **************************************************************************************************
534 : !> \brief ...
535 : !> \param param ...
536 : !> \param gfn_type ...
537 : !> \param xtb_control ...
538 : ! **************************************************************************************************
539 58 : SUBROUTINE xtb_spinpol_ext(param, gfn_type, xtb_control)
540 : TYPE(xtb_atom_type), POINTER :: param
541 : INTEGER, INTENT(IN) :: gfn_type
542 : TYPE(xtb_control_type), INTENT(IN), POINTER :: xtb_control
543 :
544 : INTEGER :: i
545 :
546 58 : SELECT CASE (gfn_type)
547 : CASE (0)
548 0 : CPABORT("gfn_type = 0: No spin polarisation possible!")
549 : CASE (1)
550 : ! OK
551 : CASE (2)
552 0 : CPABORT("gfn_type = 2 not yet supported")
553 : CASE DEFAULT
554 58 : CPABORT("Wrong gfn_type")
555 : END SELECT
556 :
557 58 : IF (param%defined) THEN
558 58 : IF (ASSOCIATED(xtb_control%spinpol_type)) THEN
559 12 : DO i = 1, SIZE(xtb_control%spinpol_type)
560 12 : IF (xtb_control%spinpol_type(i) == param%z) THEN
561 6 : param%wall(1, 1) = xtb_control%spinpol_vals(1, i)
562 6 : param%wall(1, 2) = xtb_control%spinpol_vals(2, i)
563 6 : param%wall(2, 2) = xtb_control%spinpol_vals(3, i)
564 6 : param%wall(1, 3) = xtb_control%spinpol_vals(4, i)
565 6 : param%wall(2, 3) = xtb_control%spinpol_vals(5, i)
566 6 : param%wall(3, 3) = xtb_control%spinpol_vals(6, i)
567 6 : param%wall(2, 1) = param%wall(1, 2)
568 6 : param%wall(3, 1) = param%wall(1, 3)
569 6 : param%wall(3, 2) = param%wall(2, 3)
570 6 : EXIT
571 : END IF
572 : END DO
573 : END IF
574 : END IF
575 :
576 58 : END SUBROUTINE xtb_spinpol_ext
577 :
578 : ! **************************************************************************************************
579 : !> \brief Read atom parameters for xTB Hamiltonian from input file
580 : !> \param param ...
581 : ! **************************************************************************************************
582 2256 : SUBROUTINE xtb_parameters_set(param)
583 :
584 : TYPE(xtb_atom_type), POINTER :: param
585 :
586 : INTEGER :: i, is, l, na
587 : REAL(KIND=dp), DIMENSION(5) :: kp
588 :
589 2256 : IF (param%defined) THEN
590 : ! AO to shell pointer
591 : ! AO to l-qn pointer
592 2254 : na = 0
593 7396 : DO is = 1, param%nshell
594 5142 : l = param%lval(is)
595 18558 : DO i = 1, 2*l + 1
596 11162 : na = na + 1
597 11162 : param%nao(na) = is
598 16304 : param%lao(na) = l
599 : END DO
600 : END DO
601 : !
602 2254 : i = param%z
603 : ! Electronegativity
604 2254 : param%electronegativity = eneg(i)
605 2254 : IF (param%en == 0.0_dp) param%en = eneg(i)
606 : ! covalent radius
607 2254 : param%rcov = crad(i)*bohr
608 : ! shell occupances
609 13524 : param%occupation(:) = occupation(:, i)
610 : ! check for consistency
611 13524 : IF (ABS(param%zeff - SUM(param%occupation)) > 1.E-10_dp) THEN
612 0 : CALL cp_abort(__LOCATION__, "Element <"//TRIM(param%aname)//"> has inconsistent shell occupations")
613 : END IF
614 : ! orbital energies [evolt] -> [a.u.]
615 13524 : param%hen = param%hen/evolt
616 : ! some forgotten scaling parameters (not in orig. paper)
617 2254 : param%xgamma = 0.1_dp*param%xgamma
618 13524 : param%kpoly(:) = 0.01_dp*param%kpoly(:)
619 13524 : param%kappa(:) = 0.1_dp*param%kappa(:)
620 : ! we have 1/6 g * q**3 (not 1/3)
621 2254 : param%xgamma = -2.0_dp*param%xgamma
622 : ! we need kpoly in shell order
623 13524 : kp(:) = param%kpoly(:)
624 13524 : param%kpoly(:) = 0.0_dp
625 7396 : DO is = 1, param%nshell
626 5142 : l = param%lval(is)
627 7396 : param%kpoly(is) = kp(l + 1)
628 : END DO
629 : ! kx
630 2254 : param%kx = 0.1_dp*param%kx
631 2254 : IF (param%kx < -5._dp) THEN
632 : ! use defaults
633 4462 : SELECT CASE (param%z)
634 : CASE DEFAULT
635 2220 : param%kx = 0.0_dp
636 : CASE (35) ! Br
637 12 : param%kx = 0.1_dp*0.381742_dp
638 : CASE (53) ! I
639 10 : param%kx = 0.1_dp*0.321944_dp
640 : CASE (85) ! At
641 2242 : param%kx = 0.1_dp*0.220000_dp
642 : END SELECT
643 : END IF
644 : ! chmax
645 2254 : param%chmax = clmt(i)
646 : END IF
647 :
648 2256 : END SUBROUTINE xtb_parameters_set
649 :
650 : ! **************************************************************************************************
651 : !> \brief ...
652 : !> \param param ...
653 : !> \param gto_basis_set ...
654 : !> \param ngauss ...
655 : ! **************************************************************************************************
656 2254 : SUBROUTINE init_xtb_basis(param, gto_basis_set, ngauss)
657 :
658 : TYPE(xtb_atom_type), POINTER :: param
659 : TYPE(gto_basis_set_type), POINTER :: gto_basis_set
660 : INTEGER, INTENT(IN) :: ngauss
661 :
662 2254 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: symbol
663 : INTEGER :: i, nshell
664 2254 : INTEGER, DIMENSION(:), POINTER :: lq, nq
665 2254 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet
666 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
667 :
668 2254 : IF (ASSOCIATED(param)) THEN
669 2254 : IF (param%defined) THEN
670 2254 : NULLIFY (sto_basis_set)
671 2254 : CALL allocate_sto_basis_set(sto_basis_set)
672 2254 : nshell = param%nshell
673 :
674 6762 : ALLOCATE (symbol(1:nshell))
675 7396 : symbol = ""
676 7396 : DO i = 1, nshell
677 2254 : SELECT CASE (param%lval(i))
678 : CASE (0)
679 2766 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "S"
680 : CASE (1)
681 1742 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "P"
682 : CASE (2)
683 634 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "D"
684 : CASE (3)
685 0 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "F"
686 : CASE DEFAULT
687 5142 : CPABORT('BASIS SET OUT OF RANGE (lval)')
688 : END SELECT
689 : END DO
690 :
691 2254 : IF (nshell > 0) THEN
692 13524 : ALLOCATE (nq(nshell), lq(nshell), zet(nshell))
693 14792 : nq(1:nshell) = param%nval(1:nshell)
694 14792 : lq(1:nshell) = param%lval(1:nshell)
695 14792 : zet(1:nshell) = param%zeta(1:nshell)
696 : CALL set_sto_basis_set(sto_basis_set, name=param%aname, nshell=nshell, symbol=symbol, &
697 2254 : nq=nq, lq=lq, zet=zet)
698 2254 : CALL create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss=ngauss, ortho=.TRUE.)
699 : END IF
700 :
701 : ! this will remove the allocated arrays
702 2254 : CALL deallocate_sto_basis_set(sto_basis_set)
703 2254 : DEALLOCATE (symbol, nq, lq, zet)
704 : END IF
705 :
706 : ELSE
707 0 : CPABORT("The pointer param is not associated")
708 : END IF
709 :
710 2254 : END SUBROUTINE init_xtb_basis
711 :
712 : ! **************************************************************************************************
713 : !> \brief ...
714 : !> \param za ...
715 : !> \param zb ...
716 : !> \param xtb_control ...
717 : !> \return ...
718 : ! **************************************************************************************************
719 27826 : FUNCTION xtb_set_kab(za, zb, xtb_control) RESULT(kab)
720 :
721 : INTEGER, INTENT(IN) :: za, zb
722 : TYPE(xtb_control_type), INTENT(IN), POINTER :: xtb_control
723 : REAL(KIND=dp) :: kab
724 :
725 : INTEGER :: j, z
726 : LOGICAL :: custom
727 :
728 27826 : kab = 1.0_dp
729 27826 : custom = .FALSE.
730 :
731 27826 : IF (xtb_control%kab_nval > 0) THEN
732 28 : DO j = 1, xtb_control%kab_nval
733 : IF ((za == xtb_control%kab_types(1, j) .AND. &
734 28 : zb == xtb_control%kab_types(2, j)) .OR. &
735 : (za == xtb_control%kab_types(2, j) .AND. &
736 0 : zb == xtb_control%kab_types(1, j))) THEN
737 28 : custom = .TRUE.
738 28 : kab = xtb_control%kab_vals(j)
739 : EXIT
740 : END IF
741 : END DO
742 : END IF
743 :
744 : IF (.NOT. custom) THEN
745 27798 : IF (za == 1 .OR. zb == 1) THEN
746 : ! hydrogen
747 15788 : z = za + zb - 1
748 3904 : SELECT CASE (z)
749 : CASE (1)
750 3904 : kab = 0.96_dp
751 : CASE (5)
752 0 : kab = 0.95_dp
753 : CASE (7)
754 760 : kab = 1.04_dp
755 : CASE (28)
756 0 : kab = 0.90_dp
757 : CASE (75)
758 0 : kab = 0.80_dp
759 : CASE (78)
760 15788 : kab = 0.80_dp
761 : END SELECT
762 12010 : ELSE IF (za == 5 .OR. zb == 5) THEN
763 : ! Boron
764 0 : z = za + zb - 5
765 0 : SELECT CASE (z)
766 : CASE (15)
767 0 : kab = 0.97_dp
768 : END SELECT
769 12010 : ELSE IF (za == 7 .OR. zb == 7) THEN
770 : ! Nitrogen
771 2008 : z = za + zb - 7
772 0 : SELECT CASE (z)
773 : CASE (14)
774 : !xtb orig code parameter file
775 : ! in the paper this is Kab for B-Si
776 2008 : kab = 1.01_dp
777 : END SELECT
778 10002 : ELSE IF (za > 20 .AND. za < 30) THEN
779 : ! 3d
780 178 : IF (zb > 20 .AND. zb < 30) THEN
781 : ! 3d
782 : kab = 1.10_dp
783 86 : ELSE IF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
784 : ! 4d/5d/4f
785 0 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
786 : END IF
787 9824 : ELSE IF ((za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
788 : ! 4d/5d/4f
789 30 : IF (zb > 20 .AND. zb < 30) THEN
790 : ! 3d
791 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
792 30 : ELSE IF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
793 : ! 4d/5d/4f
794 28 : kab = 1.20_dp
795 : END IF
796 : END IF
797 : END IF
798 :
799 27826 : END FUNCTION xtb_set_kab
800 :
801 : ! **************************************************************************************************
802 : !> \brief ...
803 : !> \param atag ...
804 : !> \param nshell ...
805 : !> \param nval ...
806 : !> \param lval ...
807 : !> \return ...
808 : ! **************************************************************************************************
809 2254 : SUBROUTINE xtb_get_shells(atag, nshell, nval, lval)
810 : CHARACTER(len=*) :: atag
811 : INTEGER :: nshell
812 : INTEGER, DIMENSION(:) :: nval, lval
813 :
814 : CHARACTER(LEN=1) :: ltag
815 : CHARACTER(LEN=10) :: aotag
816 : INTEGER :: i, j
817 :
818 2254 : aotag = ADJUSTL(TRIM(atag))
819 2254 : nshell = LEN(TRIM(aotag))/2
820 7396 : DO i = 1, nshell
821 5142 : j = (i - 1)*2 + 1
822 5142 : READ (aotag(j:j), FMT="(i1)") nval(i)
823 5142 : READ (aotag(j + 1:j + 1), FMT="(A1)") ltag
824 5142 : CALL uppercase(ltag)
825 2254 : SELECT CASE (ltag)
826 : CASE ("S")
827 2766 : lval(i) = 0
828 : CASE ("P")
829 1742 : lval(i) = 1
830 : CASE ("D")
831 5142 : lval(i) = 2
832 : CASE DEFAULT
833 : END SELECT
834 : END DO
835 :
836 2254 : END SUBROUTINE xtb_get_shells
837 :
838 : ! **************************************************************************************************
839 : !> \brief ...
840 : !> \param z ...
841 : !> \return ...
842 : ! **************************************************************************************************
843 10232 : FUNCTION metal(z) RESULT(ismetal)
844 : INTEGER :: z
845 : LOGICAL :: ismetal
846 :
847 10232 : SELECT CASE (z)
848 : CASE DEFAULT
849 10110 : ismetal = .TRUE.
850 : CASE (1:2, 6:10, 14:18, 32:36, 50:54, 82:86)
851 10232 : ismetal = .FALSE.
852 : END SELECT
853 :
854 10232 : END FUNCTION metal
855 :
856 : ! **************************************************************************************************
857 : !> \brief ...
858 : !> \param z ...
859 : !> \return ...
860 : ! **************************************************************************************************
861 10110 : FUNCTION early3d(z) RESULT(isearly3d)
862 : INTEGER :: z
863 : LOGICAL :: isearly3d
864 :
865 10110 : isearly3d = .FALSE.
866 10110 : IF (z >= 21 .AND. z <= 24) isearly3d = .TRUE.
867 :
868 10110 : END FUNCTION early3d
869 :
870 : ! **************************************************************************************************
871 : !> \brief ...
872 : !> \param za ...
873 : !> \param zb ...
874 : !> \return ...
875 : ! **************************************************************************************************
876 12308 : FUNCTION pp_gfn0(za, zb) RESULT(pparm)
877 : INTEGER :: za, zb
878 : REAL(KIND=dp) :: pparm
879 :
880 12308 : pparm = 1.0_dp
881 12308 : IF ((za > 20 .AND. za < 30) .OR. (za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
882 470 : IF ((zb > 20 .AND. zb < 30) .OR. (zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
883 220 : pparm = 1.1_dp
884 220 : IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
885 : IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
886 26 : pparm = 0.9_dp
887 : END IF
888 : END IF
889 : END IF
890 : END IF
891 :
892 12308 : END FUNCTION pp_gfn0
893 :
894 : END MODULE xtb_parameters
895 :
|
