1subroutine test_gfn0_sp 2 use xtb_mctc_accuracy, only : wp 3 use xtb_mctc_io, only : stdout 4 5 use assertion 6 7 use xtb_mctc_systools 8 use xtb_type_options 9 10 use xtb_type_molecule 11 use xtb_type_wavefunction 12 use xtb_type_basisset 13 use xtb_type_param 14 use xtb_type_data 15 use xtb_type_environment 16 use xtb_solv_gbsa 17 18 use xtb_setparam 19 use xtb_basis 20 use xtb_peeq 21 use xtb_readparam 22 use xtb_paramset 23 24 use xtb_xtb_data 25 use xtb_xtb_gfn0 26 27 implicit none 28 real(wp),parameter :: thr = 1.0e-7_wp 29 real(wp),parameter :: thr2 = 1.0e-5_wp 30 integer, parameter :: nat = 8 31 integer, parameter :: at(nat) = [7,15,9,9,1,1,1,1] 32 real(wp),parameter :: xyz(3,nat) = reshape(& 33 &[1.50040286526241_wp,-2.88219140061585_wp, 0.00000000000000_wp, & 34 & 0.00000000000000_wp, 1.90142164792207_wp, 0.00000000000000_wp, & 35 &-0.02649585010919_wp,-5.05651406856634_wp, 0.00000000000000_wp, & 36 &-1.39762821979929_wp, 4.65236211997835_wp, 0.00000000000000_wp, & 37 & 2.62205170116282_wp,-3.14316635901963_wp, 1.53958066744940_wp, & 38 &-1.46489869067775_wp, 0.78885483581631_wp, 1.94964934855945_wp, & 39 & 2.62205170116282_wp,-3.14316635901963_wp,-1.53958066744940_wp, & 40 &-1.46489869067775_wp, 0.78885483581631_wp,-1.94964934855945_wp],& 41 & shape(xyz)) 42 real(wp),parameter :: et = 300.0_wp 43 integer, parameter :: maxiter = 50 44 integer, parameter :: prlevel = 2 45 logical, parameter :: lgrad = .true. 46 logical, parameter :: restart = .false. 47 real(wp),parameter :: acc = 1.0_wp 48 49 type(TEnvironment) :: env 50 type(TMolecule) :: mol 51 type(scc_results) :: res 52 type(TBasisset) :: basis 53 type(TWavefunction) :: wfn 54 type(TxTBData) :: xtbData 55 type(TBorn), allocatable :: gbsa 56 57 real(wp) :: etot,egap,sigma(3,3) 58 real(wp), allocatable :: g(:,:) 59 character(len=:),allocatable :: fnv 60 integer :: ipar 61 62 type(TxTBParameter) :: globpar 63 logical :: okpar,okbas,exist,diff 64 65 call init(env) 66 67 gfn_method = 0 68 69 call init(mol, at, xyz) 70 71 wfn%nel = idint(sum(mol%z)) 72 wfn%nopen = 0 73 74 allocate( g(3,mol%n), source = 0.0_wp ) 75 76 call use_parameterset('param_gfn0-xtb.txt',globpar,xtbData,okpar) 77 !call assert(okpar) 78 call rdpath(env%xtbpath,'param_gfn0-xtb.txt',fnv,exist) 79 ! maybe the user provides a local parameter file, this was always 80 ! an option in `xtb', so we will give it a try 81 if (.not.exist) fnv = 'param_gfn0-xtb.txt' 82 call open_file(ipar,fnv,'r') 83 if (ipar.eq.-1) then 84 ! at this point there is no chance to recover from this error 85 ! THEREFORE, we have to kill the program 86 call env%error("Parameter file '"//fnv//"' not found!") 87 call terminate(1) 88 return 89 endif 90 call readParam(env,ipar,globpar,xtbData,.true.) 91 call close_file(ipar) 92 93 call newBasisset(xtbData,mol%n,mol%at,basis,okbas) 94 call assert(okbas) 95 96 call assert_eq(basis%nshell,17) 97 call assert_eq(basis%nbf, 30) 98 call assert_eq(basis%nao, 29) 99 100 call wfn%allocate(mol%n,basis%nshell,basis%nao) 101 wfn%q = mol%chrg/real(mol%n,wp) 102 103 g = 0.0_wp 104 105 call peeq(env,mol,wfn,basis,xtbData,gbsa, & 106 & egap,et,prlevel,lgrad,.false.,acc,etot,g,sigma,res) 107 108 call assert(res%converged) 109 110 call assert_close(res%e_total,-15.271927490420_wp,thr) 111 call assert_close(res%gnorm,0.55600868499186E-01_wp,thr) 112 ! value in electron volt 113 call assert_close(res%hl_gap, 4.4910991783546_wp,1.0e-4_wp) 114 call assert_close(res%e_elec,-15.238855850629_wp,thr) 115 ! es and aes are not welldefined at this accuracy level 116 call assert_close(res%e_es, -0.091088841288_wp,thr*10) 117 call assert_close(res%e_disp, -0.001331408911_wp,thr) 118 call assert_close(res%e_rep, 0.071830283799_wp,thr) 119 call assert_close(res%e_xb, -0.012481673391_wp,thr) 120 121 call mol%deallocate 122 call wfn%deallocate 123 call basis%deallocate 124 125 call terminate(afail) 126end subroutine test_gfn0_sp 127 128subroutine test_gfn0_api 129 use xtb_mctc_accuracy, only : wp 130 use xtb_mctc_io, only : stdout 131 use assertion 132 133 use xtb_type_options 134 use xtb_type_molecule 135 use xtb_type_restart 136 use xtb_type_data 137 use xtb_type_param 138 use xtb_type_environment 139 140 use xtb_pbc_tools 141 142 use xtb_xtb_calculator, only : TxTBCalculator 143 use xtb_main_setup, only : newXTBCalculator, newWavefunction 144 145 implicit none 146 147 real(wp),parameter :: thr = 1.0e-7_wp 148 integer, parameter :: nat = 7 149 integer, parameter :: at(nat) = [51,1,1,1,1,16,1] 150 real(wp),parameter :: xyz(3,nat) = reshape(& 151 [-0.33839323285714_wp,-2.52035427985714_wp, 0.00000000000000_wp, & 152 & 0.87126950614286_wp,-4.36564525285714_wp, 2.28275498100000_wp, & 153 & 0.87126950614286_wp,-4.36564525285714_wp,-2.28275498100000_wp, & 154 &-3.07540243785714_wp,-4.12867550385714_wp, 0.00000000000000_wp, & 155 & 2.16321847414286_wp, 6.15526246914286_wp, 0.00000000000000_wp, & 156 &-0.33839323285714_wp, 5.87053934414286_wp, 0.00000000000000_wp, & 157 &-0.15356858285714_wp, 3.35451847614286_wp, 0.00000000000000_wp], shape(xyz)) 158 type(peeq_options),parameter :: opt = peeq_options( & 159 & prlevel = 2, ccm = .false., acc = 1.0_wp, etemp = 300.0_wp, grad = .true. ) 160 161 type(TMolecule) :: mol 162 type(TEnvironment) :: env 163 type(TRestart) :: chk 164 type(scc_results) :: res 165 type(TxTBCalculator) :: calc 166 167 real(wp) :: energy 168 real(wp) :: hl_gap 169 real(wp) :: sigma(3,3) 170 real(wp),allocatable :: gradient(:,:) 171 172 ! setup the environment variables 173 call init(env) 174 175 call init(mol, at, xyz) 176 177 allocate(gradient(3,mol%n)) 178 energy = 0.0_wp 179 gradient = 0.0_wp 180 181 call newXTBCalculator(env, mol, calc, method=0) 182 call env%checkpoint("failed setup") 183 call newWavefunction(env, mol, calc, chk) 184 185 call calc%singlepoint(env, mol, chk, 2, .false., energy, gradient, sigma, & 186 & hl_gap, res) 187 188 call assert_close(hl_gap, 5.5384029314207_wp,thr) 189 call assert_close(energy,-8.6908532561691_wp,thr) 190 call assert_close(norm2(gradient),0.24709483673929E-01_wp,thr) 191 192 call assert_close(gradient(1,3), 0.12460501676405E-02_wp,thr) 193 call assert_close(gradient(2,1),-0.20307939298595E-01_wp,thr) 194 call assert_close(gradient(1,6),-0.29075269403059E-02_wp,thr) 195 call assert_close(gradient(3,5), 0.00000000000000E+00_wp,thr) 196 197 call terminate(afail) 198 199end subroutine test_gfn0_api 200 201subroutine test_gfn0_api_srb 202 use xtb_mctc_accuracy, only : wp 203 use xtb_mctc_io, only : stdout 204 use assertion 205 206 use xtb_type_options 207 use xtb_type_molecule 208 use xtb_type_data 209 use xtb_type_restart 210 use xtb_type_param 211 use xtb_type_environment 212 213 use xtb_pbc_tools 214 215 use xtb_xtb_calculator, only : TxTBCalculator 216 use xtb_main_setup, only : newXTBCalculator, newWavefunction 217 218 implicit none 219 220 real(wp),parameter :: thr = 1.0e-7_wp 221 integer, parameter :: nat = 24 222 integer, parameter :: at(nat) = [6,7,6,7,6,6,6,8,7,6,8,7,6,6, & 223 & 1,1,1,1,1,1,1,1,1,1] 224 real(wp),parameter :: xyz(3,nat) = reshape(& 225 &[ 2.02799738646442_wp, 0.09231312124713_wp, -0.14310895950963_wp, & 226 & 4.75011007621000_wp, 0.02373496014051_wp, -0.14324124033844_wp, & 227 & 6.33434307654413_wp, 2.07098865582721_wp, -0.14235306905930_wp, & 228 & 8.72860718071825_wp, 1.38002919517619_wp, -0.14265542523943_wp, & 229 & 8.65318821103610_wp, -1.19324866489847_wp, -0.14231527453678_wp, & 230 & 6.23857175648671_wp, -2.08353643730276_wp, -0.14218299370797_wp, & 231 & 5.63266886875962_wp, -4.69950321056008_wp, -0.13940509630299_wp, & 232 & 3.44931709749015_wp, -5.48092386085491_wp, -0.14318454855466_wp, & 233 & 7.77508917214346_wp, -6.24427872938674_wp, -0.13107140408805_wp, & 234 & 10.30229550927022_wp, -5.39739796609292_wp, -0.13672168520430_wp, & 235 & 12.07410272485492_wp, -6.91573621641911_wp, -0.13666499342053_wp, & 236 & 10.70038521493902_wp, -2.79078533715849_wp, -0.14148379504141_wp, & 237 & 13.24597858727017_wp, -1.76969072232377_wp, -0.14218299370797_wp, & 238 & 7.40891694074004_wp, -8.95905928176407_wp, -0.11636933482904_wp, & 239 & 1.38702118184179_wp, 2.05575746325296_wp, -0.14178615122154_wp, & 240 & 1.34622199478497_wp, -0.86356704498496_wp, 1.55590600570783_wp, & 241 & 1.34624089204623_wp, -0.86133716815647_wp, -1.84340893849267_wp, & 242 & 5.65596919189118_wp, 4.00172183859480_wp, -0.14131371969009_wp, & 243 & 14.67430918222276_wp, -3.26230980007732_wp, -0.14344911021228_wp, & 244 & 13.50897177220290_wp, -0.60815166181684_wp, 1.54898960808727_wp, & 245 & 13.50780014200488_wp, -0.60614855212345_wp, -1.83214617078268_wp, & 246 & 5.41408424778406_wp, -9.49239668625902_wp, -0.11022772492007_wp, & 247 & 8.31919801555568_wp, -9.74947502841788_wp, 1.56539243085954_wp, & 248 & 8.31511620712388_wp, -9.76854236502758_wp, -1.79108242206824_wp],& 249 & shape(xyz)) 250 type(peeq_options),parameter :: opt = peeq_options( & 251 & prlevel = 2, ccm = .false., acc = 1.0_wp, etemp = 300.0_wp, grad = .true. ) 252 253 type(TMolecule) :: mol 254 type(TEnvironment) :: env 255 type(TRestart) :: chk 256 type(scc_results) :: res 257 type(TxTBCalculator) :: calc 258 259 real(wp) :: energy 260 real(wp) :: hl_gap 261 real(wp) :: sigma(3,3) 262 real(wp),allocatable :: gradient(:,:) 263 264 ! setup the environment variables 265 call init(env) 266 267 call init(mol, at, xyz) 268 269 allocate(gradient(3,mol%n)) 270 energy = 0.0_wp 271 gradient = 0.0_wp 272 273 call newXTBCalculator(env, mol, calc, method=0) 274 call env%checkpoint("failed setup") 275 call newWavefunction(env, mol, calc, chk) 276 277 call calc%singlepoint(env, mol, chk, 2, .false., energy, gradient, sigma, & 278 & hl_gap, res) 279 280 call assert_close(hl_gap, 3.1192454818777_wp,thr) 281 call assert_close(energy,-40.908850360158_wp,thr) 282 call assert_close(norm2(gradient),0.83573853247451E-01_wp,thr) 283 284 call assert_close(gradient(1, 3), 0.29870640452823E-01_wp,thr) 285 call assert_close(gradient(2, 1), 0.67498165159819E-03_wp,thr) 286 call assert_close(gradient(1, 6),-0.54145817368542E-02_wp,thr) 287 call assert_close(gradient(3, 5),-0.39879097860049E-04_wp,thr) 288 call assert_close(gradient(1, 8),-0.63393710255449E-02_wp,thr) 289 call assert_close(gradient(2, 6), 0.19478157220852E-01_wp,thr) 290 call assert_close(gradient(1,20), 0.10070932785359E-03_wp,thr) 291 call assert_close(gradient(3,10),-0.54834076982200E-04_wp,thr) 292 293 call terminate(afail) 294 295end subroutine test_gfn0_api_srb 296 297 298subroutine test_gfn0_mindless_basic 299 use assertion 300 use xtb_mctc_accuracy, only : wp 301 use xtb_test_molstock, only : getMolecule 302 303 use xtb_type_molecule 304 use xtb_type_param 305 use xtb_type_pcem 306 use xtb_type_data, only : scc_results 307 use xtb_type_environment, only : TEnvironment, init 308 use xtb_type_restart, only : TRestart 309 310 use xtb_xtb_calculator, only : TxTBCalculator 311 use xtb_main_setup, only : newXTBCalculator, newWavefunction 312 313 implicit none 314 315 real(wp), parameter :: thr = 1.0e-8_wp 316 317 type(TEnvironment) :: env 318 type(TMolecule) :: mol 319 type(TRestart) :: chk 320 type(TxTBCalculator) :: calc 321 type(scc_results) :: res 322 323 integer :: iMol 324 logical :: exitRun 325 real(wp) :: energy, hl_gap, sigma(3, 3) 326 real(wp), allocatable :: gradient(:, :) 327 328 character(len=*), parameter :: mindless(10) = [& 329 & "mindless01", "mindless02", "mindless03", "mindless04", "mindless05", & 330 & "mindless06", "mindless07", "mindless08", "mindless09", "mindless10"] 331 real(wp), parameter :: ref_energies(10) = & 332 &[-28.497588460796_wp, -23.254292745575_wp, -23.060461949336_wp, & 333 & -21.822909508325_wp, -27.551996281113_wp, -18.814024099835_wp, & 334 & -31.195804700948_wp, -27.317394883750_wp, -20.258074646863_wp, & 335 & -25.051348388657_wp] 336 real(wp), parameter :: ref_gnorms(10) = & 337 &[0.05890968017315_wp, 0.04174928560148_wp, 0.05768263390495_wp, & 338 & 0.06591502692791_wp, 0.08382457849620_wp, 0.09995152858607_wp, & 339 & 0.08343577915160_wp, 0.08653790675487_wp, 0.06328938755566_wp, & 340 & 0.06269021365940_wp] 341 real(wp), parameter :: ref_hlgaps(10) = & 342 &[4.9355818428641_wp, 0.5817815942841_wp, 2.4534066058746_wp, & 343 & 0.8490082064427_wp, 1.4811842635615_wp, 0.9873933294559_wp, & 344 & 4.4138879602731_wp, 1.1871662062975_wp, 1.9816407664953_wp, & 345 & 0.3432086163775_wp] 346 347 call init(env) 348 do iMol = 1, 10 349 if (afail > 0) exit 350 351 call getMolecule(mol, mindless(iMol)) 352 353 if (allocated(gradient)) deallocate(gradient) 354 allocate(gradient(3, len(mol))) 355 356 call newXTBCalculator(env, mol, calc, method=0) 357 call newWavefunction(env, mol, calc, chk) 358 359 call env%check(exitRun) 360 call assert(.not.exitRun) 361 if (exitRun) exit 362 363 call calc%singlepoint(env, mol, chk, 2, .false., energy, gradient, sigma, & 364 & hl_gap, res) 365 366 call env%check(exitRun) 367 call assert(.not.exitRun) 368 if (exitRun) exit 369 370 call assert_close(energy, ref_energies(iMol), thr) 371 call assert_close(norm2(gradient), ref_gnorms(iMol), thr) 372 call assert_close(hl_gap, ref_hlgaps(iMol), thr) 373 374 end do 375 376 call terminate(afail) 377 378end subroutine test_gfn0_mindless_basic 379 380 381subroutine test_gfn0_mindless_solvation 382 use assertion 383 use xtb_mctc_accuracy, only : wp 384 use xtb_test_molstock, only : getMolecule 385 386 use xtb_type_molecule 387 use xtb_type_param 388 use xtb_type_pcem 389 use xtb_type_data, only : scc_results 390 use xtb_type_environment, only : TEnvironment, init 391 use xtb_type_restart, only : TRestart 392 393 use xtb_xtb_calculator, only : TxTBCalculator 394 use xtb_main_setup, only : newXTBCalculator, newWavefunction, addSolvationModel 395 use xtb_solv_input, only : TSolvInput 396 use xtb_solv_kernel, only : gbKernel 397 398 implicit none 399 400 real(wp), parameter :: thr = 1.0e-8_wp 401 402 type(TEnvironment) :: env 403 type(TMolecule) :: mol 404 type(TRestart) :: chk 405 type(TxTBCalculator) :: calc 406 type(scc_results) :: res 407 408 integer :: iMol 409 logical :: exitRun 410 real(wp) :: energy, hl_gap, sigma(3, 3) 411 real(wp), allocatable :: gradient(:, :) 412 413 character(len=*), parameter :: mindless(10) = [& 414 & "mindless01", "mindless02", "mindless03", "mindless04", "mindless05", & 415 & "mindless06", "mindless07", "mindless08", "mindless09", "mindless10"] 416 character(len=*), parameter :: solvents(10) = [character(len=20) ::& 417 & "h2o", "chcl3", "thf", "acetonitrile", "toluene", & 418 & "ch2cl2", "ether", "methanol", "cs2", "dmso"] 419 real(wp), parameter :: ref_energies(10) = & 420 &[-28.510632030412_wp, -23.266701495294_wp, -23.086125389488_wp, & 421 & -21.822303641621_wp, -27.566928060995_wp, -18.832541288587_wp, & 422 & -31.204815373848_wp, -27.323931773594_wp, -20.270903985453_wp, & 423 & -25.053546379618_wp] 424 real(wp), parameter :: ref_gnorms(10) = & 425 &[0.05788277333095_wp, 0.04117813956764_wp, 0.05838285007764_wp, & 426 & 0.06675003613700_wp, 0.08237342790744_wp, 0.10215599601359_wp, & 427 & 0.08394783095013_wp, 0.086507373560174_wp, 0.06283764829270_wp, & 428 & 0.06263571386482_wp] 429 real(wp), parameter :: ref_hlgaps(10) = & 430 &[4.9494565569973_wp, 0.5809054942177_wp, 2.4549059279536_wp, & 431 & 0.8522897920347_wp, 1.4827703858899_wp, 0.9710672494430_wp, & 432 & 4.3955269926200_wp, 1.1868023047842_wp, 1.9777105277691_wp, & 433 & 0.3654558060973_wp] 434 435 436 call init(env) 437 do iMol = 1, 10 438 if (afail > 0) exit 439 440 call getMolecule(mol, mindless(iMol)) 441 442 if (allocated(gradient)) deallocate(gradient) 443 allocate(gradient(3, len(mol))) 444 445 call newXTBCalculator(env, mol, calc, method=0) 446 call newWavefunction(env, mol, calc, chk) 447 call addSolvationModel(env, calc, TSolvInput(solvent=trim(solvents(iMol)), & 448 & alpb=mod(iMol, 2)==0, kernel=gbKernel%still)) 449 450 call env%check(exitRun) 451 call assert(.not.exitRun) 452 if (exitRun) exit 453 454 call calc%singlepoint(env, mol, chk, 2, .false., energy, gradient, sigma, & 455 & hl_gap, res) 456 457 call env%check(exitRun) 458 call assert(.not.exitRun) 459 if (exitRun) exit 460 461 call assert_close(energy, ref_energies(iMol), thr) 462 call assert_close(norm2(gradient), ref_gnorms(iMol), thr) 463 call assert_close(hl_gap, ref_hlgaps(iMol), thr) 464 465 end do 466 467 call terminate(afail) 468 469end subroutine test_gfn0_mindless_solvation 470