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