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