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_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, 2.96_dp, 3.00_dp, & ! 36
73 : 0.82_dp, 0.95_dp, 1.22_dp, 1.33_dp, 1.60_dp, 2.16_dp, 1.90_dp, 2.20_dp, &
74 : 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
75 : 0.79_dp, 0.89_dp, 1.10_dp, &
76 : 1.12_dp, 1.13_dp, 1.14_dp, 1.15_dp, 1.17_dp, 1.18_dp, 1.20_dp, 1.21_dp, &
77 : 1.22_dp, 1.23_dp, 1.24_dp, 1.25_dp, 1.26_dp, 1.27_dp, & ! Lanthanides
78 : 1.30_dp, 1.50_dp, 2.36_dp, 1.90_dp, 2.20_dp, 2.20_dp, 2.28_dp, 2.54_dp, &
79 : 2.00_dp, 2.04_dp, 2.33_dp, 2.02_dp, 2.00_dp, 2.20_dp, 2.20_dp, & ! 86
80 : 0.70_dp, 0.89_dp, 1.10_dp, &
81 : 1.30_dp, 1.50_dp, 1.38_dp, 1.36_dp, 1.28_dp, 1.30_dp, 1.30_dp, 1.30_dp, &
82 : 1.30_dp, 1.30_dp, 1.30_dp, 1.30_dp, 1.30_dp, 1.50_dp, & ! Actinides
83 : 1.50_dp, 1.50_dp, 1.50_dp/)
84 : !&>
85 :
86 : !&<
87 : ! Shell occupation
88 : INTEGER, DIMENSION(1:5, 0:nelem) :: occupation = RESHAPE((/0,0,0,0,0, & ! 0
89 : 1,0,0,0,0, 2,0,0,0,0, & ! 2
90 : 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
91 : 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
92 : 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, &
93 : 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
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, & ! 54
96 : 1,0,0,0,0, 2,0,0,0,0, 2,0,1,0,0, &
97 : 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, &
98 : 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
99 : 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, &
100 : 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)
101 : 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, & !
102 : 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, &
103 : 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
104 : 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0/), (/5, nelem+1/))
105 : !&>
106 :
107 : !&<
108 : ! COVALENT RADII
109 : ! based on "Atomic Radii of the Elements," M. Mantina, R. Valero, C. J. Cramer, and D. G. Truhlar,
110 : ! in CRC Handbook of Chemistry and Physics, 91st Edition (2010-2011),
111 : ! edited by W. M. Haynes (CRC Press, Boca Raton, FL, 2010), pages 9-49-9-50;
112 : ! corrected Nov. 17, 2010 for the 92nd edition.
113 : REAL(KIND=dp), DIMENSION(0:nelem), &
114 : PARAMETER, PRIVATE :: crad = (/0.00_dp, & ! 0
115 : 0.32_dp, 0.37_dp, & ! 2
116 : 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
117 : 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
118 : 2.00_dp, 1.74_dp, 1.59_dp, 1.48_dp, 1.44_dp, 1.30_dp, 1.29_dp, 1.24_dp, &
119 : 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
120 : 2.15_dp, 1.90_dp, 1.76_dp, 1.64_dp, 1.56_dp, 1.46_dp, 1.38_dp, 1.36_dp, &
121 : 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
122 : 2.38_dp, 2.06_dp, 1.94_dp, &
123 : 1.84_dp, 1.90_dp, 1.88_dp, 1.86_dp, 1.85_dp, 1.83_dp, 1.82_dp, 1.81_dp, &
124 : 1.80_dp, 1.79_dp, 1.77_dp, 1.77_dp, 1.78_dp, 1.74_dp, & ! Lanthanides
125 : 1.64_dp, 1.58_dp, 1.50_dp, 1.41_dp, 1.36_dp, 1.32_dp, 1.30_dp, 1.30_dp, &
126 : 1.32_dp, 1.44_dp, 1.45_dp, 1.50_dp, 1.42_dp, 1.48_dp, 1.46_dp, & ! 86
127 : 2.42_dp, 2.11_dp, 2.01_dp, &
128 : 1.90_dp, 1.84_dp, 1.83_dp, 1.80_dp, 1.80_dp, 1.51_dp, 0.96_dp, 1.54_dp, &
129 : 1.83_dp, 1.50_dp, 1.50_dp, 1.50_dp, 1.50_dp, 1.50_dp, & ! Actinides
130 : 1.50_dp, 1.50_dp, 1.50_dp/)
131 : !&>
132 :
133 : !&<
134 : ! Charge Limits (Mulliken)
135 : REAL(KIND=dp), DIMENSION(0:nelem), &
136 : PARAMETER, PRIVATE :: clmt = (/0.00_dp, & ! 0
137 : 1.05_dp, 1.25_dp, & ! 2
138 : 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
139 : 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
140 : 1.05_dp, 2.05_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
141 : 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
142 : 1.05_dp, 2.05_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
143 : 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
144 : 1.05_dp, 2.05_dp, 3.00_dp, &
145 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, &
146 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, & ! Lanthanides
147 : 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, &
148 : 2.50_dp, 2.50_dp, 3.50_dp, 3.50_dp, 3.50_dp, 1.25_dp, 1.00_dp, & ! 86
149 : 1.05_dp, 2.05_dp, 3.00_dp, &
150 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, &
151 : 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, 3.00_dp, & ! Actinides
152 : 3.00_dp, 3.00_dp, 3.00_dp/)
153 : !&>
154 :
155 : ! *** Global parameters ***
156 :
157 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_parameters'
158 :
159 : ! *** Public data types ***
160 :
161 : PUBLIC :: xtb_parameters_init, xtb_parameters_set, init_xtb_basis, xtb_set_kab
162 : PUBLIC :: metal, early3d, pp_gfn0
163 :
164 : CONTAINS
165 :
166 : ! **************************************************************************************************
167 : !> \brief ...
168 : !> \param param ...
169 : !> \param gfn_type ...
170 : !> \param element_symbol ...
171 : !> \param parameter_file_path ...
172 : !> \param parameter_file_name ...
173 : !> \param para_env ...
174 : ! **************************************************************************************************
175 2000 : SUBROUTINE xtb_parameters_init(param, gfn_type, element_symbol, &
176 : parameter_file_path, parameter_file_name, &
177 : para_env)
178 :
179 : TYPE(xtb_atom_type), POINTER :: param
180 : INTEGER, INTENT(IN) :: gfn_type
181 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
182 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
183 : TYPE(mp_para_env_type), POINTER :: para_env
184 :
185 3402 : SELECT CASE (gfn_type)
186 : CASE (0)
187 : CALL xtb0_parameters_init(param, element_symbol, parameter_file_path, &
188 1402 : parameter_file_name, para_env)
189 : CASE (1)
190 : CALL xtb1_parameters_init(param, element_symbol, parameter_file_path, &
191 598 : parameter_file_name, para_env)
192 : CASE (2)
193 0 : CPABORT("gfn_type = 2 not yet supported")
194 : CASE DEFAULT
195 2000 : CPABORT("Wrong gfn_type")
196 : END SELECT
197 :
198 2000 : END SUBROUTINE xtb_parameters_init
199 :
200 : ! **************************************************************************************************
201 : !> \brief ...
202 : !> \param param ...
203 : !> \param element_symbol ...
204 : !> \param parameter_file_path ...
205 : !> \param parameter_file_name ...
206 : !> \param para_env ...
207 : ! **************************************************************************************************
208 2804 : SUBROUTINE xtb0_parameters_init(param, element_symbol, parameter_file_path, parameter_file_name, &
209 : para_env)
210 :
211 : TYPE(xtb_atom_type), POINTER :: param
212 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
213 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
214 : TYPE(mp_para_env_type), POINTER :: para_env
215 :
216 : CHARACTER(len=2) :: esym
217 : CHARACTER(len=default_string_length) :: aname, atag, filename
218 : INTEGER :: i, l, zin, znum
219 : LOGICAL :: at_end, found
220 : TYPE(cp_parser_type) :: parser
221 :
222 1402 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(parameter_file_name))
223 1402 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
224 1402 : found = .FALSE.
225 : znum = 0
226 1402 : CALL get_ptable_info(element_symbol, znum)
227 : DO
228 : at_end = .FALSE.
229 545534 : CALL parser_get_next_line(parser, 1, at_end)
230 545534 : IF (at_end) EXIT
231 545534 : CALL parser_get_object(parser, aname)
232 545534 : CALL uppercase(aname)
233 545534 : IF (aname == "$Z") THEN
234 26836 : CALL parser_get_object(parser, zin)
235 26836 : IF (zin == znum) THEN
236 25396 : found = .TRUE.
237 : DO
238 25396 : CALL parser_get_next_line(parser, 1, at_end)
239 25396 : IF (at_end) THEN
240 0 : CPABORT("Incomplete xTB parameter file")
241 : END IF
242 25396 : CALL parser_get_object(parser, aname)
243 25396 : CALL uppercase(aname)
244 1402 : SELECT CASE (aname)
245 : CASE ("AO")
246 1402 : CALL parser_get_object(parser, atag)
247 1402 : CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
248 : CASE ("LEV")
249 4732 : DO i = 1, param%nshell
250 4732 : CALL parser_get_object(parser, param%hen(i))
251 : END DO
252 : CASE ("EXP")
253 4732 : DO i = 1, param%nshell
254 4732 : CALL parser_get_object(parser, param%zeta(i))
255 : END DO
256 : CASE ("EN")
257 578 : CALL parser_get_object(parser, param%en)
258 : CASE ("GAM")
259 1402 : CALL parser_get_object(parser, param%eta)
260 : CASE ("KQAT2")
261 1402 : CALL parser_get_object(parser, param%kqat2)
262 : CASE ("KCNS")
263 1402 : CALL parser_get_object(parser, param%kcn(1))
264 1402 : param%kcn(1) = param%kcn(1)*0.1_dp !from orig xtb code
265 : CASE ("KCNP")
266 1204 : CALL parser_get_object(parser, param%kcn(2))
267 1204 : param%kcn(2) = param%kcn(2)*0.1_dp !from orig xtb code
268 : CASE ("KCND")
269 526 : CALL parser_get_object(parser, param%kcn(3))
270 526 : param%kcn(3) = param%kcn(3)*0.1_dp !from orig xtb code
271 : CASE ("REPA")
272 1402 : CALL parser_get_object(parser, param%alpha)
273 : CASE ("REPB")
274 1402 : CALL parser_get_object(parser, param%zneff)
275 : CASE ("POLYS")
276 1402 : CALL parser_get_object(parser, param%kpoly(1))
277 : CASE ("POLYP")
278 1204 : CALL parser_get_object(parser, param%kpoly(2))
279 : CASE ("POLYD")
280 526 : CALL parser_get_object(parser, param%kpoly(3))
281 : CASE ("KQS")
282 1402 : CALL parser_get_object(parser, param%kq(1))
283 : CASE ("KQP")
284 1204 : CALL parser_get_object(parser, param%kq(2))
285 : CASE ("KQD")
286 526 : CALL parser_get_object(parser, param%kq(3))
287 : CASE ("XI")
288 1402 : CALL parser_get_object(parser, param%xi)
289 : CASE ("KAPPA")
290 1402 : CALL parser_get_object(parser, param%kappa0)
291 : CASE ("ALPG")
292 1402 : CALL parser_get_object(parser, param%alpg)
293 : CASE ("$END")
294 0 : EXIT
295 : CASE DEFAULT
296 25396 : CPABORT("Unknown parameter in xTB file")
297 : END SELECT
298 : END DO
299 : ELSE
300 : CYCLE
301 : END IF
302 : EXIT
303 : END IF
304 : END DO
305 1402 : IF (found) THEN
306 1402 : param%typ = "STANDARD"
307 1402 : param%symbol = element_symbol
308 1402 : param%defined = .TRUE.
309 1402 : param%z = znum
310 1402 : param%aname = ptable(znum)%name
311 4732 : param%lmax = MAXVAL(param%lval(1:param%nshell))
312 1402 : param%natorb = 0
313 4732 : DO i = 1, param%nshell
314 3330 : l = param%lval(i)
315 4732 : param%natorb = param%natorb + (2*l + 1)
316 : END DO
317 1402 : param%zeff = zval(znum)
318 : ELSE
319 0 : esym = element_symbol
320 0 : CALL uppercase(esym)
321 0 : IF ("X " == esym) THEN
322 0 : param%typ = "GHOST"
323 0 : param%symbol = element_symbol
324 0 : param%defined = .FALSE.
325 0 : param%z = 0
326 0 : param%aname = "X "
327 0 : param%lmax = 0
328 0 : param%natorb = 0
329 0 : param%nshell = 0
330 0 : param%zeff = 0.0_dp
331 : ELSE
332 0 : param%defined = .FALSE.
333 : CALL cp_warn(__LOCATION__, "xTB parameters for element "//element_symbol// &
334 0 : " were not found in the parameter file "//ADJUSTL(TRIM(filename)))
335 : END IF
336 : END IF
337 1402 : CALL parser_release(parser)
338 :
339 4206 : END SUBROUTINE xtb0_parameters_init
340 :
341 : ! **************************************************************************************************
342 : !> \brief ...
343 : !> \param param ...
344 : !> \param element_symbol ...
345 : !> \param parameter_file_path ...
346 : !> \param parameter_file_name ...
347 : !> \param para_env ...
348 : ! **************************************************************************************************
349 1196 : SUBROUTINE xtb1_parameters_init(param, element_symbol, parameter_file_path, parameter_file_name, &
350 : para_env)
351 :
352 : TYPE(xtb_atom_type), POINTER :: param
353 : CHARACTER(LEN=2), INTENT(IN) :: element_symbol
354 : CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, parameter_file_name
355 : TYPE(mp_para_env_type), POINTER :: para_env
356 :
357 : CHARACTER(len=2) :: esym
358 : CHARACTER(len=default_string_length) :: aname, atag, filename
359 : INTEGER :: i, l, zin, znum
360 : LOGICAL :: at_end, found
361 : TYPE(cp_parser_type) :: parser
362 :
363 598 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(parameter_file_name))
364 598 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
365 598 : found = .FALSE.
366 : znum = 0
367 598 : CALL get_ptable_info(element_symbol, znum)
368 : DO
369 : at_end = .FALSE.
370 72152 : CALL parser_get_next_line(parser, 1, at_end)
371 72152 : IF (at_end) EXIT
372 72150 : CALL parser_get_object(parser, aname)
373 72150 : CALL uppercase(aname)
374 72152 : IF (aname == "$Z") THEN
375 4784 : CALL parser_get_object(parser, zin)
376 4784 : IF (zin == znum) THEN
377 5770 : found = .TRUE.
378 : DO
379 5770 : CALL parser_get_next_line(parser, 1, at_end)
380 5770 : IF (at_end) THEN
381 0 : CPABORT("Incomplete xTB parameter file")
382 : END IF
383 5770 : CALL parser_get_object(parser, aname)
384 5770 : CALL uppercase(aname)
385 596 : SELECT CASE (aname)
386 : CASE ("AO")
387 596 : CALL parser_get_object(parser, atag)
388 596 : CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
389 : CASE ("LEV")
390 1862 : DO i = 1, param%nshell
391 1862 : CALL parser_get_object(parser, param%hen(i))
392 : END DO
393 : CASE ("EXP")
394 1862 : DO i = 1, param%nshell
395 1862 : CALL parser_get_object(parser, param%zeta(i))
396 : END DO
397 : CASE ("GAM")
398 596 : CALL parser_get_object(parser, param%eta)
399 : CASE ("GAM3")
400 368 : CALL parser_get_object(parser, param%xgamma)
401 : CASE ("CXB")
402 12 : CALL parser_get_object(parser, param%kx)
403 : CASE ("REPA")
404 596 : CALL parser_get_object(parser, param%alpha)
405 : CASE ("REPB")
406 596 : CALL parser_get_object(parser, param%zneff)
407 : CASE ("POLYS")
408 368 : CALL parser_get_object(parser, param%kpoly(1))
409 : CASE ("POLYP")
410 368 : CALL parser_get_object(parser, param%kpoly(2))
411 : CASE ("POLYD")
412 74 : CALL parser_get_object(parser, param%kpoly(3))
413 : CASE ("LPARP")
414 360 : CALL parser_get_object(parser, param%kappa(2))
415 : CASE ("LPARD")
416 48 : CALL parser_get_object(parser, param%kappa(3))
417 : CASE ("$END")
418 0 : EXIT
419 : CASE DEFAULT
420 5770 : CPABORT("Unknown parameter in xTB file")
421 : END SELECT
422 : END DO
423 : ELSE
424 : CYCLE
425 : END IF
426 : EXIT
427 : END IF
428 : END DO
429 598 : IF (found) THEN
430 596 : param%typ = "STANDARD"
431 596 : param%symbol = element_symbol
432 596 : param%defined = .TRUE.
433 596 : param%z = znum
434 596 : param%aname = ptable(znum)%name
435 1862 : param%lmax = MAXVAL(param%lval(1:param%nshell))
436 596 : param%natorb = 0
437 1862 : DO i = 1, param%nshell
438 1266 : l = param%lval(i)
439 1862 : param%natorb = param%natorb + (2*l + 1)
440 : END DO
441 596 : param%zeff = zval(znum)
442 : ELSE
443 2 : esym = element_symbol
444 2 : CALL uppercase(esym)
445 2 : IF ("X " == esym) THEN
446 2 : param%typ = "GHOST"
447 2 : param%symbol = element_symbol
448 2 : param%defined = .FALSE.
449 2 : param%z = 0
450 2 : param%aname = "X "
451 2 : param%lmax = 0
452 2 : param%natorb = 0
453 2 : param%nshell = 0
454 2 : param%zeff = 0.0_dp
455 : ELSE
456 0 : param%defined = .FALSE.
457 : CALL cp_warn(__LOCATION__, "xTB parameters for element "//element_symbol// &
458 0 : " were not found in the parameter file "//ADJUSTL(TRIM(filename)))
459 : END IF
460 : END IF
461 598 : CALL parser_release(parser)
462 :
463 1794 : END SUBROUTINE xtb1_parameters_init
464 :
465 : ! **************************************************************************************************
466 : !> \brief Read atom parameters for xTB Hamiltonian from input file
467 : !> \param param ...
468 : ! **************************************************************************************************
469 2000 : SUBROUTINE xtb_parameters_set(param)
470 :
471 : TYPE(xtb_atom_type), POINTER :: param
472 :
473 : INTEGER :: i, is, l, na
474 : REAL(KIND=dp), DIMENSION(5) :: kp
475 :
476 2000 : IF (param%defined) THEN
477 : ! AO to shell pointer
478 : ! AO to l-qn pointer
479 1998 : na = 0
480 6594 : DO is = 1, param%nshell
481 4596 : l = param%lval(is)
482 16734 : DO i = 1, 2*l + 1
483 10140 : na = na + 1
484 10140 : param%nao(na) = is
485 14736 : param%lao(na) = l
486 : END DO
487 : END DO
488 : !
489 1998 : i = param%z
490 : ! Electronegativity
491 1998 : param%electronegativity = eneg(i)
492 1998 : IF (param%en == 0.0_dp) param%en = eneg(i)
493 : ! covalent radius
494 1998 : param%rcov = crad(i)*bohr
495 : ! shell occupances
496 11988 : param%occupation(:) = occupation(:, i)
497 : ! check for consistency
498 11988 : IF (ABS(param%zeff - SUM(param%occupation)) > 1.E-10_dp) THEN
499 0 : CALL cp_abort(__LOCATION__, "Element <"//TRIM(param%aname)//"> has inconsistent shell occupations")
500 : END IF
501 : ! orbital energies [evolt] -> [a.u.]
502 11988 : param%hen = param%hen/evolt
503 : ! some forgotten scaling parameters (not in orig. paper)
504 1998 : param%xgamma = 0.1_dp*param%xgamma
505 11988 : param%kpoly(:) = 0.01_dp*param%kpoly(:)
506 11988 : param%kappa(:) = 0.1_dp*param%kappa(:)
507 : ! we have 1/6 g * q**3 (not 1/3)
508 1998 : param%xgamma = -2.0_dp*param%xgamma
509 : ! we need kpoly in shell order
510 11988 : kp(:) = param%kpoly(:)
511 11988 : param%kpoly(:) = 0.0_dp
512 6594 : DO is = 1, param%nshell
513 4596 : l = param%lval(is)
514 6594 : param%kpoly(is) = kp(l + 1)
515 : END DO
516 : ! kx
517 1998 : param%kx = 0.1_dp*param%kx
518 1998 : IF (param%kx < -5._dp) THEN
519 : ! use defaults
520 3950 : SELECT CASE (param%z)
521 : CASE DEFAULT
522 1964 : param%kx = 0.0_dp
523 : CASE (35) ! Br
524 12 : param%kx = 0.1_dp*0.381742_dp
525 : CASE (53) ! I
526 10 : param%kx = 0.1_dp*0.321944_dp
527 : CASE (85) ! At
528 1986 : param%kx = 0.1_dp*0.220000_dp
529 : END SELECT
530 : END IF
531 : ! chmax
532 1998 : param%chmax = clmt(i)
533 : END IF
534 :
535 2000 : END SUBROUTINE xtb_parameters_set
536 :
537 : ! **************************************************************************************************
538 : !> \brief ...
539 : !> \param param ...
540 : !> \param gto_basis_set ...
541 : !> \param ngauss ...
542 : ! **************************************************************************************************
543 1998 : SUBROUTINE init_xtb_basis(param, gto_basis_set, ngauss)
544 :
545 : TYPE(xtb_atom_type), POINTER :: param
546 : TYPE(gto_basis_set_type), POINTER :: gto_basis_set
547 : INTEGER, INTENT(IN) :: ngauss
548 :
549 1998 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: symbol
550 : INTEGER :: i, nshell
551 1998 : INTEGER, DIMENSION(:), POINTER :: lq, nq
552 1998 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet
553 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
554 :
555 1998 : IF (ASSOCIATED(param)) THEN
556 1998 : IF (param%defined) THEN
557 1998 : NULLIFY (sto_basis_set)
558 1998 : CALL allocate_sto_basis_set(sto_basis_set)
559 1998 : nshell = param%nshell
560 :
561 5994 : ALLOCATE (symbol(1:nshell))
562 6594 : symbol = ""
563 6594 : DO i = 1, nshell
564 1998 : SELECT CASE (param%lval(i))
565 : CASE (0)
566 2424 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "S"
567 : CASE (1)
568 1572 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "P"
569 : CASE (2)
570 600 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "D"
571 : CASE (3)
572 0 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "F"
573 : CASE DEFAULT
574 4596 : CPABORT('BASIS SET OUT OF RANGE (lval)')
575 : END SELECT
576 : END DO
577 :
578 1998 : IF (nshell > 0) THEN
579 11988 : ALLOCATE (nq(nshell), lq(nshell), zet(nshell))
580 13188 : nq(1:nshell) = param%nval(1:nshell)
581 13188 : lq(1:nshell) = param%lval(1:nshell)
582 13188 : zet(1:nshell) = param%zeta(1:nshell)
583 : CALL set_sto_basis_set(sto_basis_set, name=param%aname, nshell=nshell, symbol=symbol, &
584 1998 : nq=nq, lq=lq, zet=zet)
585 1998 : CALL create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss=ngauss, ortho=.TRUE.)
586 : END IF
587 :
588 : ! this will remove the allocated arrays
589 1998 : CALL deallocate_sto_basis_set(sto_basis_set)
590 1998 : DEALLOCATE (symbol, nq, lq, zet)
591 : END IF
592 :
593 : ELSE
594 0 : CPABORT("The pointer param is not associated")
595 : END IF
596 :
597 1998 : END SUBROUTINE init_xtb_basis
598 :
599 : ! **************************************************************************************************
600 : !> \brief ...
601 : !> \param za ...
602 : !> \param zb ...
603 : !> \param xtb_control ...
604 : !> \return ...
605 : ! **************************************************************************************************
606 21834 : FUNCTION xtb_set_kab(za, zb, xtb_control) RESULT(kab)
607 :
608 : INTEGER, INTENT(IN) :: za, zb
609 : TYPE(xtb_control_type), INTENT(IN), POINTER :: xtb_control
610 : REAL(KIND=dp) :: kab
611 :
612 : INTEGER :: j, z
613 : LOGICAL :: custom
614 :
615 21834 : kab = 1.0_dp
616 21834 : custom = .FALSE.
617 :
618 21834 : IF (xtb_control%kab_nval .GT. 0) THEN
619 28 : DO j = 1, xtb_control%kab_nval
620 : IF ((za == xtb_control%kab_types(1, j) .AND. &
621 28 : zb == xtb_control%kab_types(2, j)) .OR. &
622 : (za == xtb_control%kab_types(2, j) .AND. &
623 0 : zb == xtb_control%kab_types(1, j))) THEN
624 28 : custom = .TRUE.
625 28 : kab = xtb_control%kab_vals(j)
626 28 : EXIT
627 : END IF
628 : END DO
629 : END IF
630 :
631 21834 : IF (.NOT. custom) THEN
632 21806 : IF (za == 1 .OR. zb == 1) THEN
633 : ! hydrogen
634 12482 : z = za + zb - 1
635 3102 : SELECT CASE (z)
636 : CASE (1)
637 3102 : kab = 0.96_dp
638 : CASE (5)
639 0 : kab = 0.95_dp
640 : CASE (7)
641 760 : kab = 1.04_dp
642 : CASE (28)
643 0 : kab = 0.90_dp
644 : CASE (75)
645 0 : kab = 0.80_dp
646 : CASE (78)
647 12482 : kab = 0.80_dp
648 : END SELECT
649 9324 : ELSEIF (za == 5 .OR. zb == 5) THEN
650 : ! Boron
651 0 : z = za + zb - 5
652 0 : SELECT CASE (z)
653 : CASE (15)
654 0 : kab = 0.97_dp
655 : END SELECT
656 9324 : ELSEIF (za == 7 .OR. zb == 7) THEN
657 : ! Nitrogen
658 2008 : z = za + zb - 7
659 0 : SELECT CASE (z)
660 : CASE (14)
661 : !xtb orig code parameter file
662 : ! in the paper this is Kab for B-Si
663 2008 : kab = 1.01_dp
664 : END SELECT
665 7316 : ELSEIF (za > 20 .AND. za < 30) THEN
666 : ! 3d
667 96 : IF (zb > 20 .AND. zb < 30) THEN
668 : ! 3d
669 : kab = 1.10_dp
670 76 : ELSEIF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
671 : ! 4d/5d/4f
672 0 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
673 : END IF
674 7220 : ELSEIF ((za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
675 : ! 4d/5d/4f
676 30 : IF (zb > 20 .AND. zb < 30) THEN
677 : ! 3d
678 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
679 30 : ELSEIF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
680 : ! 4d/5d/4f
681 28 : kab = 1.20_dp
682 : END IF
683 : END IF
684 : END IF
685 :
686 21834 : END FUNCTION xtb_set_kab
687 :
688 : ! **************************************************************************************************
689 : !> \brief ...
690 : !> \param atag ...
691 : !> \param nshell ...
692 : !> \param nval ...
693 : !> \param lval ...
694 : !> \return ...
695 : ! **************************************************************************************************
696 1998 : SUBROUTINE xtb_get_shells(atag, nshell, nval, lval)
697 : CHARACTER(len=*) :: atag
698 : INTEGER :: nshell
699 : INTEGER, DIMENSION(:) :: nval, lval
700 :
701 : CHARACTER(LEN=1) :: ltag
702 : CHARACTER(LEN=10) :: aotag
703 : INTEGER :: i, j
704 :
705 1998 : aotag = ADJUSTL(TRIM(atag))
706 1998 : nshell = LEN(TRIM(aotag))/2
707 6594 : DO i = 1, nshell
708 4596 : j = (i - 1)*2 + 1
709 4596 : READ (aotag(j:j), FMT="(i1)") nval(i)
710 4596 : READ (aotag(j + 1:j + 1), FMT="(A1)") ltag
711 4596 : CALL uppercase(ltag)
712 1998 : SELECT CASE (ltag)
713 : CASE ("S")
714 2424 : lval(i) = 0
715 : CASE ("P")
716 1572 : lval(i) = 1
717 : CASE ("D")
718 4596 : lval(i) = 2
719 : CASE DEFAULT
720 : END SELECT
721 : END DO
722 :
723 1998 : END SUBROUTINE xtb_get_shells
724 :
725 : ! **************************************************************************************************
726 : !> \brief ...
727 : !> \param z ...
728 : !> \return ...
729 : ! **************************************************************************************************
730 7872 : FUNCTION metal(z) RESULT(ismetal)
731 : INTEGER :: z
732 : LOGICAL :: ismetal
733 :
734 7872 : SELECT CASE (z)
735 : CASE DEFAULT
736 7822 : ismetal = .TRUE.
737 : CASE (1:2, 6:10, 14:18, 32:36, 50:54, 82:86)
738 7872 : ismetal = .FALSE.
739 : END SELECT
740 :
741 7872 : END FUNCTION metal
742 :
743 : ! **************************************************************************************************
744 : !> \brief ...
745 : !> \param z ...
746 : !> \return ...
747 : ! **************************************************************************************************
748 7822 : FUNCTION early3d(z) RESULT(isearly3d)
749 : INTEGER :: z
750 : LOGICAL :: isearly3d
751 :
752 7822 : isearly3d = .FALSE.
753 7822 : IF (z >= 21 .AND. z <= 24) isearly3d = .TRUE.
754 :
755 7822 : END FUNCTION early3d
756 :
757 : ! **************************************************************************************************
758 : !> \brief ...
759 : !> \param za ...
760 : !> \param zb ...
761 : !> \return ...
762 : ! **************************************************************************************************
763 8742 : FUNCTION pp_gfn0(za, zb) RESULT(pparm)
764 : INTEGER :: za, zb
765 : REAL(KIND=dp) :: pparm
766 :
767 8742 : pparm = 1.0_dp
768 8742 : IF ((za > 20 .AND. za < 30) .OR. (za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
769 470 : IF ((zb > 20 .AND. zb < 30) .OR. (zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
770 220 : pparm = 1.1_dp
771 220 : IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
772 : IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
773 26 : pparm = 0.9_dp
774 : END IF
775 : END IF
776 : END IF
777 : END IF
778 :
779 8742 : END FUNCTION pp_gfn0
780 :
781 : END MODULE xtb_parameters
782 :
|
