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 : !&<
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 :
189 : PUBLIC :: xtb_parameters_init, xtb_parameters_set, init_xtb_basis, xtb_set_kab
190 : PUBLIC :: xtb_spinpol_init, xtb_spinpol_ext
191 : PUBLIC :: metal, early3d, pp_gfn0
192 :
193 : CONTAINS
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 2296 : 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 3770 : SELECT CASE (gfn_type)
215 : CASE (0)
216 : CALL xtb0_parameters_init(param, element_symbol, parameter_file_path, &
217 1474 : parameter_file_name, para_env)
218 : CASE (1)
219 : CALL xtb1_parameters_init(param, element_symbol, parameter_file_path, &
220 822 : parameter_file_name, para_env)
221 : CASE (2)
222 0 : CPABORT("gfn_type = 2 not yet supported")
223 : CASE DEFAULT
224 2296 : CPABORT("Wrong gfn_type")
225 : END SELECT
226 :
227 2296 : 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 2948 : 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 1474 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(parameter_file_name))
252 1474 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
253 1474 : found = .FALSE.
254 : znum = 0
255 1474 : CALL get_ptable_info(element_symbol, znum)
256 : DO
257 : at_end = .FALSE.
258 554742 : CALL parser_get_next_line(parser, 1, at_end)
259 554742 : IF (at_end) EXIT
260 554742 : CALL parser_get_object(parser, aname)
261 554742 : CALL uppercase(aname)
262 554742 : IF (aname == "$Z") THEN
263 27298 : CALL parser_get_object(parser, zin)
264 27298 : IF (zin == znum) THEN
265 26628 : found = .TRUE.
266 : DO
267 26628 : CALL parser_get_next_line(parser, 1, at_end)
268 26628 : IF (at_end) THEN
269 0 : CPABORT("Incomplete xTB parameter file")
270 : END IF
271 26628 : CALL parser_get_object(parser, aname)
272 26628 : CALL uppercase(aname)
273 1474 : SELECT CASE (aname)
274 : CASE ("AO")
275 1474 : CALL parser_get_object(parser, atag)
276 1474 : CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
277 : CASE ("LEV")
278 4960 : DO i = 1, param%nshell
279 4960 : CALL parser_get_object(parser, param%hen(i))
280 : END DO
281 : CASE ("EXP")
282 4960 : DO i = 1, param%nshell
283 4960 : CALL parser_get_object(parser, param%zeta(i))
284 : END DO
285 : CASE ("EN")
286 616 : CALL parser_get_object(parser, param%en)
287 : CASE ("GAM")
288 1474 : CALL parser_get_object(parser, param%eta)
289 : CASE ("KQAT2")
290 1474 : CALL parser_get_object(parser, param%kqat2)
291 : CASE ("KCNS")
292 1474 : CALL parser_get_object(parser, param%kcn(1))
293 1474 : param%kcn(1) = param%kcn(1)*0.1_dp !from orig xtb code
294 : CASE ("KCNP")
295 1254 : CALL parser_get_object(parser, param%kcn(2))
296 1254 : param%kcn(2) = param%kcn(2)*0.1_dp !from orig xtb code
297 : CASE ("KCND")
298 538 : CALL parser_get_object(parser, param%kcn(3))
299 538 : param%kcn(3) = param%kcn(3)*0.1_dp !from orig xtb code
300 : CASE ("REPA")
301 1474 : CALL parser_get_object(parser, param%alpha)
302 : CASE ("REPB")
303 1474 : CALL parser_get_object(parser, param%zneff)
304 : CASE ("POLYS")
305 1474 : CALL parser_get_object(parser, param%kpoly(1))
306 : CASE ("POLYP")
307 1254 : CALL parser_get_object(parser, param%kpoly(2))
308 : CASE ("POLYD")
309 538 : CALL parser_get_object(parser, param%kpoly(3))
310 : CASE ("KQS")
311 1474 : CALL parser_get_object(parser, param%kq(1))
312 : CASE ("KQP")
313 1254 : CALL parser_get_object(parser, param%kq(2))
314 : CASE ("KQD")
315 538 : CALL parser_get_object(parser, param%kq(3))
316 : CASE ("XI")
317 1474 : CALL parser_get_object(parser, param%xi)
318 : CASE ("KAPPA")
319 1474 : CALL parser_get_object(parser, param%kappa0)
320 : CASE ("ALPG")
321 1474 : CALL parser_get_object(parser, param%alpg)
322 : CASE ("$END")
323 0 : EXIT
324 : CASE DEFAULT
325 26628 : 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 1474 : IF (found) THEN
335 1474 : param%typ = "STANDARD"
336 1474 : param%symbol = element_symbol
337 1474 : param%defined = .TRUE.
338 1474 : param%z = znum
339 1474 : param%aname = ptable(znum)%name
340 4960 : param%lmax = MAXVAL(param%lval(1:param%nshell))
341 1474 : param%natorb = 0
342 4960 : DO i = 1, param%nshell
343 3486 : l = param%lval(i)
344 4960 : param%natorb = param%natorb + (2*l + 1)
345 : END DO
346 1474 : param%zeff = zval(znum)
347 : ELSE
348 0 : esym = element_symbol
349 0 : CALL uppercase(esym)
350 0 : IF ("X " == esym) THEN
351 0 : param%typ = "GHOST"
352 0 : param%symbol = element_symbol
353 0 : param%defined = .FALSE.
354 0 : param%z = 0
355 0 : param%aname = "X "
356 0 : param%lmax = 0
357 0 : param%natorb = 0
358 0 : param%nshell = 0
359 0 : param%zeff = 0.0_dp
360 : ELSE
361 0 : param%defined = .FALSE.
362 : CALL cp_warn(__LOCATION__, "xTB parameters for element "//element_symbol// &
363 0 : " were not found in the parameter file "//ADJUSTL(TRIM(filename)))
364 : END IF
365 : END IF
366 1474 : CALL parser_release(parser)
367 :
368 4422 : 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 1644 : 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 822 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(parameter_file_name))
393 822 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
394 822 : found = .FALSE.
395 : znum = 0
396 822 : CALL get_ptable_info(element_symbol, znum)
397 : DO
398 : at_end = .FALSE.
399 98646 : CALL parser_get_next_line(parser, 1, at_end)
400 98646 : IF (at_end) EXIT
401 98644 : CALL parser_get_object(parser, aname)
402 98644 : CALL uppercase(aname)
403 98646 : IF (aname == "$Z") THEN
404 6552 : CALL parser_get_object(parser, zin)
405 6552 : IF (zin == znum) THEN
406 7980 : found = .TRUE.
407 : DO
408 7980 : CALL parser_get_next_line(parser, 1, at_end)
409 7980 : IF (at_end) THEN
410 0 : CPABORT("Incomplete xTB parameter file")
411 : END IF
412 7980 : CALL parser_get_object(parser, aname)
413 7980 : CALL uppercase(aname)
414 820 : SELECT CASE (aname)
415 : CASE ("AO")
416 820 : CALL parser_get_object(parser, atag)
417 820 : CALL xtb_get_shells(atag, param%nshell, param%nval, param%lval)
418 : CASE ("LEV")
419 2566 : DO i = 1, param%nshell
420 2566 : CALL parser_get_object(parser, param%hen(i))
421 : END DO
422 : CASE ("EXP")
423 2566 : DO i = 1, param%nshell
424 2566 : CALL parser_get_object(parser, param%zeta(i))
425 : END DO
426 : CASE ("GAM")
427 820 : CALL parser_get_object(parser, param%eta)
428 : CASE ("GAM3")
429 518 : CALL parser_get_object(parser, param%xgamma)
430 : CASE ("CXB")
431 12 : CALL parser_get_object(parser, param%kx)
432 : CASE ("REPA")
433 820 : CALL parser_get_object(parser, param%alpha)
434 : CASE ("REPB")
435 820 : CALL parser_get_object(parser, param%zneff)
436 : CASE ("POLYS")
437 518 : CALL parser_get_object(parser, param%kpoly(1))
438 : CASE ("POLYP")
439 518 : CALL parser_get_object(parser, param%kpoly(2))
440 : CASE ("POLYD")
441 106 : CALL parser_get_object(parser, param%kpoly(3))
442 : CASE ("LPARP")
443 508 : CALL parser_get_object(parser, param%kappa(2))
444 : CASE ("LPARD")
445 60 : CALL parser_get_object(parser, param%kappa(3))
446 : CASE ("$END")
447 0 : EXIT
448 : CASE DEFAULT
449 7980 : 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 822 : IF (found) THEN
459 820 : param%typ = "STANDARD"
460 820 : param%symbol = element_symbol
461 820 : param%defined = .TRUE.
462 820 : param%z = znum
463 820 : param%aname = ptable(znum)%name
464 2566 : param%lmax = MAXVAL(param%lval(1:param%nshell))
465 820 : param%natorb = 0
466 2566 : DO i = 1, param%nshell
467 1746 : l = param%lval(i)
468 2566 : param%natorb = param%natorb + (2*l + 1)
469 : END DO
470 820 : param%zeff = zval(znum)
471 : ELSE
472 2 : esym = element_symbol
473 2 : CALL uppercase(esym)
474 2 : IF ("X " == esym) THEN
475 2 : param%typ = "GHOST"
476 2 : param%symbol = element_symbol
477 2 : param%defined = .FALSE.
478 2 : param%z = 0
479 2 : param%aname = "X "
480 2 : param%lmax = 0
481 2 : param%natorb = 0
482 2 : param%nshell = 0
483 2 : param%zeff = 0.0_dp
484 : ELSE
485 0 : param%defined = .FALSE.
486 : CALL cp_warn(__LOCATION__, "xTB parameters for element "//element_symbol// &
487 0 : " were not found in the parameter file "//ADJUSTL(TRIM(filename)))
488 : END IF
489 : END IF
490 822 : CALL parser_release(parser)
491 :
492 2466 : 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 116 : 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 58 : SELECT CASE (gfn_type)
519 : CASE (0)
520 0 : CPABORT("gfn_type = 0: No spin polarisation possible!")
521 : CASE (1)
522 : ! OK
523 : CASE (2)
524 0 : CPABORT("gfn_type = 2 not yet supported")
525 : CASE DEFAULT
526 0 : CPABORT("Wrong gfn_type")
527 : END SELECT
528 :
529 58 : filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(spinpol_param_file_name))
530 58 : CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env)
531 : znum = 0
532 754 : param%wall = 0.0_dp
533 58 : CALL get_ptable_info(element_symbol, znum)
534 4904 : DO
535 : at_end = .FALSE.
536 4962 : CALL parser_get_next_line(parser, 1, at_end)
537 4962 : IF (at_end) EXIT
538 4904 : CALL parser_get_object(parser, zin)
539 4962 : IF (zin == znum) THEN
540 58 : CALL parser_get_object(parser, param%wall(1, 1))
541 58 : CALL parser_get_object(parser, param%wall(1, 2))
542 58 : CALL parser_get_object(parser, param%wall(2, 2))
543 58 : CALL parser_get_object(parser, param%wall(1, 3))
544 58 : CALL parser_get_object(parser, param%wall(2, 3))
545 58 : CALL parser_get_object(parser, param%wall(3, 3))
546 58 : param%wall(2, 1) = param%wall(1, 2)
547 58 : param%wall(3, 1) = param%wall(1, 3)
548 58 : param%wall(3, 2) = param%wall(2, 3)
549 : END IF
550 : END DO
551 58 : CALL parser_release(parser)
552 :
553 174 : END SUBROUTINE xtb_spinpol_init
554 :
555 : ! **************************************************************************************************
556 : !> \brief ...
557 : !> \param param ...
558 : !> \param gfn_type ...
559 : !> \param xtb_control ...
560 : ! **************************************************************************************************
561 58 : 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 58 : SELECT CASE (gfn_type)
569 : CASE (0)
570 0 : CPABORT("gfn_type = 0: No spin polarisation possible!")
571 : CASE (1)
572 : ! OK
573 : CASE (2)
574 0 : CPABORT("gfn_type = 2 not yet supported")
575 : CASE DEFAULT
576 58 : CPABORT("Wrong gfn_type")
577 : END SELECT
578 :
579 58 : IF (param%defined) THEN
580 58 : IF (ASSOCIATED(xtb_control%spinpol_type)) THEN
581 12 : DO i = 1, SIZE(xtb_control%spinpol_type)
582 12 : IF (xtb_control%spinpol_type(i) == param%z) THEN
583 6 : param%wall(1, 1) = xtb_control%spinpol_vals(1, i)
584 6 : param%wall(1, 2) = xtb_control%spinpol_vals(2, i)
585 6 : param%wall(2, 2) = xtb_control%spinpol_vals(3, i)
586 6 : param%wall(1, 3) = xtb_control%spinpol_vals(4, i)
587 6 : param%wall(2, 3) = xtb_control%spinpol_vals(5, i)
588 6 : param%wall(3, 3) = xtb_control%spinpol_vals(6, i)
589 6 : param%wall(2, 1) = param%wall(1, 2)
590 6 : param%wall(3, 1) = param%wall(1, 3)
591 6 : param%wall(3, 2) = param%wall(2, 3)
592 6 : EXIT
593 : END IF
594 : END DO
595 : END IF
596 : END IF
597 :
598 58 : END SUBROUTINE xtb_spinpol_ext
599 :
600 : ! **************************************************************************************************
601 : !> \brief Read atom parameters for xTB Hamiltonian from input file
602 : !> \param param ...
603 : ! **************************************************************************************************
604 2296 : 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 2296 : IF (param%defined) THEN
612 : ! AO to shell pointer
613 : ! AO to l-qn pointer
614 2294 : na = 0
615 7526 : DO is = 1, param%nshell
616 5232 : l = param%lval(is)
617 18878 : DO i = 1, 2*l + 1
618 11352 : na = na + 1
619 11352 : param%nao(na) = is
620 16584 : param%lao(na) = l
621 : END DO
622 : END DO
623 : !
624 2294 : i = param%z
625 : ! Electronegativity
626 2294 : param%electronegativity = eneg(i)
627 2294 : IF (param%en == 0.0_dp) param%en = eneg(i)
628 : ! covalent radius
629 2294 : param%rcov = crad(i)*bohr
630 : ! shell occupances
631 13764 : param%occupation(:) = occupation(:, i)
632 : ! number of primitive Gaussians per shell
633 13764 : param%ngauss = 0
634 9176 : param%ngauss(1:3) = number_of_primitives(:, i)
635 : ! check for consistency
636 13764 : IF (ABS(param%zeff - SUM(param%occupation)) > 1.E-10_dp) THEN
637 0 : CALL cp_abort(__LOCATION__, "Element <"//TRIM(param%aname)//"> has inconsistent shell occupations")
638 : END IF
639 : ! orbital energies [evolt] -> [a.u.]
640 13764 : param%hen = param%hen/evolt
641 : ! some forgotten scaling parameters (not in orig. paper)
642 2294 : param%xgamma = 0.1_dp*param%xgamma
643 13764 : param%kpoly(:) = 0.01_dp*param%kpoly(:)
644 13764 : param%kappa(:) = 0.1_dp*param%kappa(:)
645 : ! we have 1/6 g * q**3 (not 1/3)
646 2294 : param%xgamma = -2.0_dp*param%xgamma
647 : ! we need kpoly in shell order
648 13764 : kp(:) = param%kpoly(:)
649 13764 : param%kpoly(:) = 0.0_dp
650 7526 : DO is = 1, param%nshell
651 5232 : l = param%lval(is)
652 7526 : param%kpoly(is) = kp(l + 1)
653 : END DO
654 : ! kx
655 2294 : param%kx = 0.1_dp*param%kx
656 2294 : IF (param%kx < -5._dp) THEN
657 : ! use defaults
658 4542 : SELECT CASE (param%z)
659 : CASE DEFAULT
660 2260 : param%kx = 0.0_dp
661 : CASE (35) ! Br
662 12 : param%kx = 0.1_dp*0.381742_dp
663 : CASE (53) ! I
664 10 : param%kx = 0.1_dp*0.321944_dp
665 : CASE (85) ! At
666 2282 : param%kx = 0.1_dp*0.220000_dp
667 : END SELECT
668 : END IF
669 : ! chmax
670 2294 : param%chmax = clmt(i)
671 : END IF
672 :
673 2296 : 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 2294 : 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 2294 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: symbol
690 : INTEGER :: i, nshell
691 2294 : INTEGER, DIMENSION(:), POINTER :: lq, nq
692 2294 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet
693 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
694 :
695 2294 : IF (ASSOCIATED(param)) THEN
696 2294 : IF (param%defined) THEN
697 2294 : NULLIFY (sto_basis_set)
698 2294 : CALL allocate_sto_basis_set(sto_basis_set)
699 2294 : nshell = param%nshell
700 :
701 6882 : ALLOCATE (symbol(1:nshell))
702 7526 : symbol = ""
703 7526 : DO i = 1, nshell
704 2294 : SELECT CASE (param%lval(i))
705 : CASE (0)
706 2816 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "S"
707 : CASE (1)
708 1772 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "P"
709 : CASE (2)
710 644 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "D"
711 : CASE (3)
712 0 : WRITE (symbol(i), '(I1,A1)') param%nval(i), "F"
713 : CASE DEFAULT
714 5232 : CPABORT('BASIS SET OUT OF RANGE (lval)')
715 : END SELECT
716 : END DO
717 :
718 2294 : IF (nshell > 0) THEN
719 13764 : ALLOCATE (nq(nshell), lq(nshell), zet(nshell))
720 15052 : nq(1:nshell) = param%nval(1:nshell)
721 15052 : lq(1:nshell) = param%lval(1:nshell)
722 15052 : zet(1:nshell) = param%zeta(1:nshell)
723 : CALL set_sto_basis_set(sto_basis_set, name=param%aname, nshell=nshell, symbol=symbol, &
724 2294 : nq=nq, lq=lq, zet=zet)
725 2294 : IF (PRESENT(ngaussflex)) THEN
726 : CALL create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss=ngauss, ortho=.TRUE., &
727 8 : ngaussflex=ngaussflex)
728 : ELSE
729 2286 : 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 2294 : CALL deallocate_sto_basis_set(sto_basis_set)
735 2294 : DEALLOCATE (symbol, nq, lq, zet)
736 : END IF
737 :
738 : ELSE
739 0 : CPABORT("The pointer param is not associated")
740 : END IF
741 :
742 2294 : END SUBROUTINE init_xtb_basis
743 :
744 : ! **************************************************************************************************
745 : !> \brief ...
746 : !> \param za ...
747 : !> \param zb ...
748 : !> \param xtb_control ...
749 : !> \return ...
750 : ! **************************************************************************************************
751 28078 : 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 28078 : kab = 1.0_dp
761 28078 : custom = .FALSE.
762 :
763 28078 : IF (xtb_control%kab_nval > 0) THEN
764 28 : DO j = 1, xtb_control%kab_nval
765 : IF ((za == xtb_control%kab_types(1, j) .AND. &
766 28 : zb == xtb_control%kab_types(2, j)) .OR. &
767 : (za == xtb_control%kab_types(2, j) .AND. &
768 0 : zb == xtb_control%kab_types(1, j))) THEN
769 28 : custom = .TRUE.
770 28 : kab = xtb_control%kab_vals(j)
771 : EXIT
772 : END IF
773 : END DO
774 : END IF
775 :
776 : IF (.NOT. custom) THEN
777 28050 : IF (za == 1 .OR. zb == 1) THEN
778 : ! hydrogen
779 15856 : z = za + zb - 1
780 3924 : SELECT CASE (z)
781 : CASE (1)
782 3924 : kab = 0.96_dp
783 : CASE (5)
784 0 : kab = 0.95_dp
785 : CASE (7)
786 760 : kab = 1.04_dp
787 : CASE (28)
788 0 : kab = 0.90_dp
789 : CASE (75)
790 0 : kab = 0.80_dp
791 : CASE (78)
792 15856 : kab = 0.80_dp
793 : END SELECT
794 12194 : ELSE IF (za == 5 .OR. zb == 5) THEN
795 : ! Boron
796 0 : z = za + zb - 5
797 0 : SELECT CASE (z)
798 : CASE (15)
799 0 : kab = 0.97_dp
800 : END SELECT
801 12194 : ELSE IF (za == 7 .OR. zb == 7) THEN
802 : ! Nitrogen
803 2008 : z = za + zb - 7
804 0 : SELECT CASE (z)
805 : CASE (14)
806 : !xtb orig code parameter file
807 : ! in the paper this is Kab for B-Si
808 2008 : kab = 1.01_dp
809 : END SELECT
810 10186 : ELSE IF (za > 20 .AND. za < 30) THEN
811 : ! 3d
812 178 : IF (zb > 20 .AND. zb < 30) THEN
813 : ! 3d
814 : kab = 1.10_dp
815 86 : ELSE IF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
816 : ! 4d/5d/4f
817 0 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
818 : END IF
819 10008 : ELSE IF ((za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
820 : ! 4d/5d/4f
821 30 : IF (zb > 20 .AND. zb < 30) THEN
822 : ! 3d
823 : kab = 0.50_dp*(1.20_dp + 1.10_dp)
824 30 : ELSE IF ((zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
825 : ! 4d/5d/4f
826 28 : kab = 1.20_dp
827 : END IF
828 : END IF
829 : END IF
830 :
831 28078 : 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 2294 : 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 2294 : aotag = ADJUSTL(TRIM(atag))
851 2294 : nshell = LEN(TRIM(aotag))/2
852 7526 : DO i = 1, nshell
853 5232 : j = (i - 1)*2 + 1
854 5232 : READ (aotag(j:j), FMT="(i1)") nval(i)
855 5232 : READ (aotag(j + 1:j + 1), FMT="(A1)") ltag
856 5232 : CALL uppercase(ltag)
857 2294 : SELECT CASE (ltag)
858 : CASE ("S")
859 2816 : lval(i) = 0
860 : CASE ("P")
861 1772 : lval(i) = 1
862 : CASE ("D")
863 5232 : lval(i) = 2
864 : CASE DEFAULT
865 : END SELECT
866 : END DO
867 :
868 2294 : END SUBROUTINE xtb_get_shells
869 :
870 : ! **************************************************************************************************
871 : !> \brief ...
872 : !> \param z ...
873 : !> \return ...
874 : ! **************************************************************************************************
875 10356 : FUNCTION metal(z) RESULT(ismetal)
876 : INTEGER :: z
877 : LOGICAL :: ismetal
878 :
879 10356 : SELECT CASE (z)
880 : CASE DEFAULT
881 10202 : ismetal = .TRUE.
882 : CASE (1:2, 6:10, 14:18, 32:36, 50:54, 82:86)
883 10356 : ismetal = .FALSE.
884 : END SELECT
885 :
886 10356 : END FUNCTION metal
887 :
888 : ! **************************************************************************************************
889 : !> \brief ...
890 : !> \param z ...
891 : !> \return ...
892 : ! **************************************************************************************************
893 10202 : FUNCTION early3d(z) RESULT(isearly3d)
894 : INTEGER :: z
895 : LOGICAL :: isearly3d
896 :
897 10202 : isearly3d = .FALSE.
898 10202 : IF (z >= 21 .AND. z <= 24) isearly3d = .TRUE.
899 :
900 10202 : END FUNCTION early3d
901 :
902 : ! **************************************************************************************************
903 : !> \brief ...
904 : !> \param za ...
905 : !> \param zb ...
906 : !> \return ...
907 : ! **************************************************************************************************
908 12308 : FUNCTION pp_gfn0(za, zb) RESULT(pparm)
909 : INTEGER :: za, zb
910 : REAL(KIND=dp) :: pparm
911 :
912 12308 : pparm = 1.0_dp
913 12308 : IF ((za > 20 .AND. za < 30) .OR. (za > 38 .AND. za < 48) .OR. (za > 56 .AND. za < 80)) THEN
914 470 : IF ((zb > 20 .AND. zb < 30) .OR. (zb > 38 .AND. zb < 48) .OR. (zb > 56 .AND. zb < 80)) THEN
915 220 : pparm = 1.1_dp
916 220 : IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
917 : IF (za == 29 .OR. za == 47 .OR. za == 79) THEN
918 26 : pparm = 0.9_dp
919 : END IF
920 : END IF
921 : END IF
922 : END IF
923 :
924 12308 : END FUNCTION pp_gfn0
925 :
926 : END MODULE xtb_parameters
|
