(git:ccc2433)
qs_dispersion_pairpot.F
Go to the documentation of this file.
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 Calculation of dispersion using pair potentials
10 !> \author JGH
11 ! **************************************************************************************************
12 MODULE qs_dispersion_pairpot
13 
14  USE atomic_kind_types, ONLY: atomic_kind_type,&
15  get_atomic_kind,&
16  get_atomic_kind_set
17  USE atprop_types, ONLY: atprop_array_init,&
18  atprop_type
19  USE bibliography, ONLY: goerigk2017,&
20  cite_reference,&
21  grimme2006,&
22  grimme2010,&
23  grimme2011
24  USE cell_types, ONLY: cell_type,&
25  get_cell,&
26  pbc,&
27  plane_distance
28  USE cp_files, ONLY: close_file,&
29  open_file
30  USE cp_log_handling, ONLY: cp_get_default_logger,&
31  cp_logger_type
32  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
33  cp_print_key_unit_nr
34  USE cp_parser_methods, ONLY: parser_get_next_line,&
35  parser_get_object
36  USE cp_parser_types, ONLY: cp_parser_type,&
37  parser_create,&
38  parser_release
39  USE input_constants, ONLY: vdw_pairpot_dftd2,&
40  vdw_pairpot_dftd3,&
41  vdw_pairpot_dftd3bj,&
42  xc_vdw_fun_none,&
43  xc_vdw_fun_pairpot
44  USE input_section_types, ONLY: section_vals_type,&
45  section_vals_val_get
46  USE kinds, ONLY: default_string_length,&
47  dp
48  USE mathconstants, ONLY: twopi
49  USE memory_utilities, ONLY: reallocate
50  USE message_passing, ONLY: mp_para_env_type
51  USE molecule_types, ONLY: molecule_type
52  USE particle_types, ONLY: particle_type
53  USE physcon, ONLY: bohr,&
54  kcalmol,&
55  kjmol
56  USE qs_dispersion_types, ONLY: dftd2_pp,&
57  dftd3_pp,&
58  qs_atom_dispersion_type,&
59  qs_dispersion_type
60  USE qs_environment_types, ONLY: get_qs_env,&
61  qs_environment_type
62  USE qs_force_types, ONLY: qs_force_type
63  USE qs_kind_types, ONLY: get_qs_kind,&
64  qs_kind_type,&
65  set_qs_kind
66  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
67  neighbor_list_iterate,&
68  neighbor_list_iterator_create,&
69  neighbor_list_iterator_p_type,&
70  neighbor_list_iterator_release,&
71  neighbor_list_set_p_type
72  USE string_utilities, ONLY: uppercase
73  USE virial_methods, ONLY: virial_pair_force
74  USE virial_types, ONLY: virial_type
75 
76 !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
77 !$ omp_get_thread_num, &
78 !$ omp_lock_kind, &
79 !$ omp_init_lock, omp_set_lock, &
80 !$ omp_unset_lock, omp_destroy_lock
81 
82 #include "./base/base_uses.f90"
83 
84  IMPLICIT NONE
85 
86  PRIVATE
87 
88  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_pairpot'
89 
90  PUBLIC :: qs_dispersion_pairpot_init, calculate_dispersion_pairpot, &
91  qs_scaling_init, qs_scaling_dftd3, qs_scaling_dftd3bj, &
92  qs_dispersion_setcn
93 
94  TYPE dcnum_type
95  INTEGER :: neighbors
96  INTEGER, DIMENSION(:), POINTER :: nlist
97  REAL(KIND=dp), DIMENSION(:), POINTER :: dvals
98  REAL(KIND=dp), DIMENSION(:, :), POINTER :: rik
99  END TYPE dcnum_type
100 
101  PUBLIC :: d3_cnumber, dcnum_type, dcnum_distribute
102 
103 ! **************************************************************************************************
104 
105 CONTAINS
106 
107 ! **************************************************************************************************
108 !> \brief ...
109 !> \param atomic_kind_set ...
110 !> \param qs_kind_set ...
111 !> \param dispersion_env ...
112 !> \param pp_section ...
113 !> \param para_env ...
114 ! **************************************************************************************************
115  SUBROUTINE qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
116  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
117  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
118  TYPE(qs_dispersion_type), POINTER :: dispersion_env
119  TYPE(section_vals_type), OPTIONAL, POINTER :: pp_section
120  TYPE(mp_para_env_type), POINTER :: para_env
121 
122  CHARACTER(len=*), PARAMETER :: routinen = 'qs_dispersion_pairpot_init'
123 
124  CHARACTER(LEN=2) :: symbol
125  CHARACTER(LEN=default_string_length) :: aname, filename
126  CHARACTER(LEN=default_string_length), &
127  DIMENSION(:), POINTER :: tmpstringlist
128  INTEGER :: elem, handle, i, ikind, j, max_elem, &
129  maxc, n_rep, nkind, nl, vdw_pp_type, &
130  vdw_type
131  INTEGER, DIMENSION(:), POINTER :: exlist
132  LOGICAL :: at_end, explicit, found, is_available
133  REAL(kind=dp) :: dum
134  TYPE(qs_atom_dispersion_type), POINTER :: disp
135 
136  CALL timeset(routinen, handle)
137 
138  vdw_type = dispersion_env%type
139  SELECT CASE (vdw_type)
140  CASE DEFAULT
141  ! do nothing
142  CASE (xc_vdw_fun_pairpot)
143  ! setup information on pair potentials
144  vdw_pp_type = dispersion_env%type
145  SELECT CASE (dispersion_env%pp_type)
146  CASE DEFAULT
147  ! do nothing
148  CASE (vdw_pairpot_dftd3, vdw_pairpot_dftd3bj)
149  !DFT-D3 Method initial setup
150  CALL cite_reference(grimme2010)
151  CALL cite_reference(grimme2011)
152  CALL cite_reference(goerigk2017)
153  max_elem = 94
154  maxc = 5
155  dispersion_env%max_elem = max_elem
156  dispersion_env%maxc = maxc
157  ALLOCATE (dispersion_env%maxci(max_elem))
158  ALLOCATE (dispersion_env%c6ab(max_elem, max_elem, maxc, maxc, 3))
159  ALLOCATE (dispersion_env%r0ab(max_elem, max_elem))
160  ALLOCATE (dispersion_env%rcov(max_elem))
161  ALLOCATE (dispersion_env%r2r4(max_elem))
162  ALLOCATE (dispersion_env%cn(max_elem))
163 
164  ! get filename of parameter file
165  filename = dispersion_env%parameter_file_name
166  CALL dftd3_c6_param(dispersion_env%c6ab, dispersion_env%maxci, filename, para_env)
167  CALL setr0ab(dispersion_env%r0ab, dispersion_env%rcov, dispersion_env%r2r4)
168  ! the default coordination numbers
169  CALL setcn(dispersion_env%cn)
170  ! scale r4/r2 values of the atoms by sqrt(Z)
171  ! sqrt is also globally close to optimum
172  ! together with the factor 1/2 this yield reasonable
173  ! c8 for he, ne and ar. for larger Z, C8 becomes too large
174  ! which effectively mimics higher R^n terms neglected due
175  ! to stability reasons
176  DO i = 1, max_elem
177  dum = 0.5_dp*dispersion_env%r2r4(i)*real(i, dp)**0.5_dp
178  ! store it as sqrt because the geom. av. is taken
179  dispersion_env%r2r4(i) = sqrt(dum)
180  END DO
181  ! parameters
182  dispersion_env%k1 = 16.0_dp
183  dispersion_env%k2 = 4._dp/3._dp
184  ! reasonable choices are between 3 and 5
185  ! this gives smoth curves with maxima around the integer values
186  ! k3=3 give for CN=0 a slightly smaller value than computed
187  ! for the free atom. This also yields to larger CN for atoms
188  ! in larger molecules but with the same chem. environment
189  ! which is physically not right
190  ! values >5 might lead to bumps in the potential
191  dispersion_env%k3 = -4._dp
192  dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*bohr
193  ! alpha default parameter
194  dispersion_env%alp = 14._dp
195  END SELECT
196  END SELECT
197 
198  nkind = SIZE(atomic_kind_set)
199  SELECT CASE (vdw_type)
200  CASE DEFAULT
201  cpabort("")
202  CASE (xc_vdw_fun_none)
203  ! we should never reach this point
204  cpabort("xc_vdw_fun_none in init routine")
205  CASE (xc_vdw_fun_pairpot)
206  ! setup information on pair potentials
207  SELECT CASE (dispersion_env%pp_type)
208  CASE DEFAULT
209  cpabort("")
210  CASE (vdw_pairpot_dftd2)
211  CALL cite_reference(grimme2006)
212  DO ikind = 1, nkind
213  CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
214  ALLOCATE (disp)
215  disp%type = dftd2_pp
216  ! get filename of parameter file
217  filename = dispersion_env%parameter_file_name
218  ! check for local parameters
219  found = .false.
220  IF (PRESENT(pp_section)) THEN
221  CALL section_vals_val_get(pp_section, "ATOMPARM", n_rep_val=n_rep)
222  DO i = 1, n_rep
223  CALL section_vals_val_get(pp_section, "ATOMPARM", i_rep_val=i, &
224  c_vals=tmpstringlist)
225  IF (trim(tmpstringlist(1)) == trim(symbol)) THEN
226  ! we assume the parameters are in atomic units!
227  READ (tmpstringlist(2), *) disp%c6
228  READ (tmpstringlist(3), *) disp%vdw_radii
229  found = .true.
230  EXIT
231  END IF
232  END DO
233  END IF
234  IF (.NOT. found) THEN
235  ! check for internal parameters
236  CALL dftd2_param(elem, disp%c6, disp%vdw_radii, found)
237  END IF
238  IF (.NOT. found) THEN
239  ! check on file
240  INQUIRE (file=filename, exist=is_available)
241  IF (is_available) THEN
242  block
243  TYPE(cp_parser_type) :: parser
244  CALL parser_create(parser, filename, para_env=para_env)
245  DO
246  at_end = .false.
247  CALL parser_get_next_line(parser, 1, at_end)
248  IF (at_end) EXIT
249  CALL parser_get_object(parser, aname)
250  IF (trim(aname) == trim(symbol)) THEN
251  CALL parser_get_object(parser, disp%c6)
252  ! we have to change the units J*nm^6*mol^-1 -> Hartree*Bohr^6
253  disp%c6 = disp%c6*1000._dp*bohr**6/kjmol
254  CALL parser_get_object(parser, disp%vdw_radii)
255  disp%vdw_radii = disp%vdw_radii*bohr
256  found = .true.
257  EXIT
258  END IF
259  END DO
260  CALL parser_release(parser)
261  END block
262  END IF
263  END IF
264  IF (found) THEN
265  disp%defined = .true.
266  ELSE
267  disp%defined = .false.
268  END IF
269  ! Check if the parameter is defined
270  IF (.NOT. disp%defined) &
271  CALL cp_abort(__location__, &
272  "Dispersion parameters for element ("//trim(symbol)//") are not defined! "// &
273  "Please provide a valid set of parameters through the input section or "// &
274  "through an external file! ")
275  CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
276  END DO
277  CASE (vdw_pairpot_dftd3, vdw_pairpot_dftd3bj)
278  DO ikind = 1, nkind
279  CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
280  ALLOCATE (disp)
281  disp%type = dftd3_pp
282  IF (elem <= 94) THEN
283  disp%defined = .true.
284  ELSE
285  disp%defined = .false.
286  END IF
287  IF (.NOT. disp%defined) &
288  CALL cp_abort(__location__, &
289  "Dispersion parameters for element ("//trim(symbol)//") are not defined! "// &
290  "Please provide a valid set of parameters through the input section or "// &
291  "through an external file! ")
292  CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
293  END DO
294 
295  IF (PRESENT(pp_section)) THEN
296  ! Check for coordination numbers
297  CALL section_vals_val_get(pp_section, "KIND_COORDINATION_NUMBERS", n_rep_val=n_rep)
298  IF (n_rep > 0) THEN
299  ALLOCATE (dispersion_env%cnkind(n_rep))
300  DO i = 1, n_rep
301  CALL section_vals_val_get(pp_section, "KIND_COORDINATION_NUMBERS", i_rep_val=i, &
302  c_vals=tmpstringlist)
303  READ (tmpstringlist(1), *) dispersion_env%cnkind(i)%cnum
304  READ (tmpstringlist(2), *) dispersion_env%cnkind(i)%kind
305  END DO
306  END IF
307  CALL section_vals_val_get(pp_section, "ATOM_COORDINATION_NUMBERS", n_rep_val=n_rep)
308  IF (n_rep > 0) THEN
309  ALLOCATE (dispersion_env%cnlist(n_rep))
310  DO i = 1, n_rep
311  CALL section_vals_val_get(pp_section, "ATOM_COORDINATION_NUMBERS", i_rep_val=i, &
312  c_vals=tmpstringlist)
313  nl = SIZE(tmpstringlist)
314  ALLOCATE (dispersion_env%cnlist(i)%atom(nl - 1))
315  dispersion_env%cnlist(i)%natom = nl - 1
316  READ (tmpstringlist(1), *) dispersion_env%cnlist(i)%cnum
317  DO j = 1, nl - 1
318  READ (tmpstringlist(j + 1), *) dispersion_env%cnlist(i)%atom(j)
319  END DO
320  END DO
321  END IF
322  ! Check for exclusion lists
323  CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND", explicit=explicit)
324  IF (explicit) THEN
325  CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND", i_vals=exlist)
326  DO j = 1, SIZE(exlist)
327  ikind = exlist(j)
328  CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp)
329  disp%defined = .false.
330  END DO
331  END IF
332  CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", n_rep_val=n_rep)
333  dispersion_env%nd3_exclude_pair = n_rep
334  IF (n_rep > 0) THEN
335  ALLOCATE (dispersion_env%d3_exclude_pair(n_rep, 2))
336  DO i = 1, n_rep
337  CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", i_rep_val=i, &
338  i_vals=exlist)
339  dispersion_env%d3_exclude_pair(i, :) = exlist
340  END DO
341  END IF
342  END IF
343 
344  END SELECT
345  END SELECT
346 
347  CALL timestop(handle)
348 
349  END SUBROUTINE qs_dispersion_pairpot_init
350 
351 ! **************************************************************************************************
352 !> \brief ...
353 !> \param atomic_kind_set ...
354 !> \param qs_kind_set ...
355 !> \param dispersion_env ...
356 !> \param para_env ...
357 ! **************************************************************************************************
358  SUBROUTINE qs_dispersion_setcn(atomic_kind_set, qs_kind_set, dispersion_env, para_env)
359  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
360  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
361  TYPE(qs_dispersion_type), POINTER :: dispersion_env
362  TYPE(mp_para_env_type), POINTER :: para_env
363 
364  CHARACTER(LEN=default_string_length) :: filename
365  INTEGER :: i, ikind, max_elem, maxc, nkind
366  REAL(kind=dp) :: dum
367  TYPE(qs_atom_dispersion_type), POINTER :: disp
368 
369  nkind = SIZE(atomic_kind_set)
370  IF (dispersion_env%type /= xc_vdw_fun_pairpot) THEN
371  DO ikind = 1, nkind
372  ALLOCATE (disp)
373  disp%type = dftd3_pp
374  disp%defined = .true.
375  CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
376  END DO
377  END IF
378 
379  max_elem = 94
380  maxc = 5
381  dispersion_env%max_elem = max_elem
382  dispersion_env%maxc = maxc
383  ALLOCATE (dispersion_env%maxci(max_elem))
384  ALLOCATE (dispersion_env%c6ab(max_elem, max_elem, maxc, maxc, 3))
385  ALLOCATE (dispersion_env%r0ab(max_elem, max_elem))
386  ALLOCATE (dispersion_env%rcov(max_elem))
387  ALLOCATE (dispersion_env%r2r4(max_elem))
388  ALLOCATE (dispersion_env%cn(max_elem))
389  ! get filename of parameter file
390  filename = dispersion_env%parameter_file_name
391  CALL dftd3_c6_param(dispersion_env%c6ab, dispersion_env%maxci, filename, para_env)
392  CALL setr0ab(dispersion_env%r0ab, dispersion_env%rcov, dispersion_env%r2r4)
393  ! the default coordination numbers
394  CALL setcn(dispersion_env%cn)
395  DO i = 1, max_elem
396  dum = 0.5_dp*dispersion_env%r2r4(i)*real(i, dp)**0.5_dp
397  dispersion_env%r2r4(i) = sqrt(dum)
398  END DO
399  dispersion_env%k1 = 16.0_dp
400  dispersion_env%k2 = 4._dp/3._dp
401  dispersion_env%k3 = -4._dp
402  dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*bohr
403  dispersion_env%alp = 14._dp
404 
405  END SUBROUTINE qs_dispersion_setcn
406 
407 ! **************************************************************************************************
408 !> \brief ...
409 !> \param scaling ...
410 !> \param vdw_section ...
411 ! **************************************************************************************************
412  SUBROUTINE qs_scaling_init(scaling, vdw_section)
413  REAL(kind=dp), INTENT(inout) :: scaling
414  TYPE(section_vals_type), POINTER :: vdw_section
415 
416  CHARACTER(LEN=default_string_length) :: functional
417 
418  CALL section_vals_val_get(vdw_section, "PAIR_POTENTIAL%REFERENCE_FUNCTIONAL", c_val=functional)
419 
420  SELECT CASE (trim(functional))
421  CASE DEFAULT
422  ! unknown functional
423  cpabort("No DFT-D2 s6 value available for this functional:"//trim(functional))
424  CASE ("BLYP")
425  scaling = 1.20_dp
426  CASE ("B3LYP")
427  scaling = 1.05_dp
428  CASE ("TPSS")
429  scaling = 1.00_dp
430  CASE ("PBE")
431  scaling = 0.75_dp
432  CASE ("PBE0")
433  scaling = 0.6_dp
434  CASE ("B2PLYP")
435  scaling = 0.55_dp
436  CASE ("BP86")
437  scaling = 1.05_dp
438  CASE ("B97")
439  scaling = 1.25_dp
440  END SELECT
441 
442  END SUBROUTINE qs_scaling_init
443 
444 ! **************************************************************************************************
445 
446 ! **************************************************************************************************
447 !> \brief ...
448 !> \param s6 ...
449 !> \param sr6 ...
450 !> \param s8 ...
451 !> \param vdw_section ...
452 ! **************************************************************************************************
453  SUBROUTINE qs_scaling_dftd3(s6, sr6, s8, vdw_section)
454 
455  REAL(kind=dp), INTENT(inout) :: s6, sr6, s8
456  TYPE(section_vals_type), POINTER :: vdw_section
457 
458  CHARACTER(LEN=default_string_length) :: functional
459 
460  CALL section_vals_val_get(vdw_section, "PAIR_POTENTIAL%REFERENCE_FUNCTIONAL", c_val=functional)
461  CALL uppercase(functional)
462  ! values for different functionals from:
463  ! https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3
464  ! L. Goerigk et al. PCCP 2017, 32147-32744, SI
465  SELECT CASE (trim(functional))
466  CASE DEFAULT
467  ! unknown functional
468  cpabort("No DFT-D3 values available for this functional:"//trim(functional))
469  CASE ("B1B95")
470  s6 = 1.000_dp
471  sr6 = 1.613_dp
472  s8 = 1.868_dp
473  CASE ("B2GPPLYP")
474  ! L. Goerigk and S. Grimme
475  ! J. Chem. Theory Comput. 2011, 7, 291-309; doi:10.1021/ct100466k
476  s6 = 0.56_dp
477  sr6 = 1.586_dp
478  s8 = 0.760_dp
479  CASE ("B2PLYP")
480  ! L. Goerigk and S. Grimme
481  ! J. Chem. Theory Comput. 2011, 7, 291-309; doi:10.1021/ct100466k
482  s6 = 0.64_dp
483  sr6 = 1.427_dp
484  s8 = 1.022_dp
485  CASE ("DSD-BLYP")
486  ! L. Goerigk and S. Grimme
487  ! J. Chem. Theory Comput. 2011, 7, 291-309; doi:10.1021/ct100466k
488  s6 = 0.50_dp
489  sr6 = 1.569_dp
490  s8 = 0.705_dp
491  CASE ("B3LYP")
492  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
493  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
494  s6 = 1.000_dp
495  sr6 = 1.261_dp
496  s8 = 1.703_dp
497  CASE ("B97-D")
498  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
499  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
500  s6 = 1.000_dp
501  sr6 = 0.892_dp
502  s8 = 0.909_dp
503  CASE ("BHLYP")
504  s6 = 1.000_dp
505  sr6 = 1.370_dp
506  s8 = 1.442_dp
507  CASE ("BLYP")
508  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
509  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
510  s6 = 1.000_dp
511  sr6 = 1.094_dp
512  s8 = 1.682_dp
513  CASE ("BP86")
514  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
515  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
516  s6 = 1.000_dp
517  sr6 = 1.139_dp
518  s8 = 1.683_dp
519  CASE ("BPBE")
520  s6 = 1.000_dp
521  sr6 = 1.087_dp
522  s8 = 2.033_dp
523  CASE ("MPWLYP")
524  s6 = 1.000_dp
525  sr6 = 1.239_dp
526  s8 = 1.098_dp
527  CASE ("PBE")
528  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
529  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
530  s6 = 1.000_dp
531  sr6 = 1.217_dp
532  s8 = 0.722_dp
533  CASE ("PBEHPBE")
534  s6 = 1.000_dp
535  sr6 = 1.5703_dp
536  s8 = 1.4010_dp
537  CASE ("PBE0")
538  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
539  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
540  s6 = 1.000_dp
541  sr6 = 1.287_dp
542  s8 = 0.928_dp
543  CASE ("PW6B95")
544  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
545  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
546  s6 = 1.000_dp
547  sr6 = 1.532_dp
548  s8 = 0.862_dp
549  CASE ("PWB6K")
550  s6 = 1.000_dp
551  sr6 = 1.660_dp
552  s8 = 0.550_dp
553  CASE ("REVPBE")
554  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
555  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
556  s6 = 1.000_dp
557  sr6 = 0.923_dp
558  s8 = 1.010_dp
559  CASE ("RPBE")
560  s6 = 1.000_dp
561  sr6 = 0.872_dp
562  s8 = 0.514_dp
563  CASE ("TPSS")
564  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
565  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
566  s6 = 1.000_dp
567  sr6 = 1.166_dp
568  s8 = 1.105_dp
569  CASE ("TPSS0")
570  ! S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
571  ! J. Chem. Phys. 132,154104 (2010); doi:10.1063/1.3382344
572  s6 = 1.000_dp
573  sr6 = 1.252_dp
574  s8 = 1.242_dp
575  CASE ("TPSSH")
576  s6 = 1.000_dp
577  sr6 = 1.223_dp
578  s8 = 1.219_dp
579  CASE ("B1LYP")
580  s6 = 1.000_dp
581  sr6 = 1.3725_dp
582  s8 = 1.9467_dp
583  CASE ("B1P86")
584  s6 = 1.000_dp
585  sr6 = 1.1815_dp
586  s8 = 1.1209_dp
587  CASE ("B3P86")
588  s6 = 1.000_dp
589  sr6 = 1.1897_dp
590  s8 = 1.1961_dp
591  CASE ("B3PW91")
592  s6 = 1.000_dp
593  sr6 = 1.176_dp
594  s8 = 1.775_dp
595  CASE ("BMK")
596  s6 = 1.000_dp
597  sr6 = 1.931_dp
598  s8 = 2.168_dp
599  CASE ("CAMB3LYP")
600  s6 = 1.000_dp
601  sr6 = 1.378_dp
602  s8 = 1.217_dp
603  CASE ("LCWPBE")
604  s6 = 1.000_dp
605  sr6 = 1.355_dp
606  s8 = 1.279_dp
607  CASE ("M052X")
608  s6 = 1.000_dp
609  sr6 = 1.417_dp
610  s8 = 0.000_dp
611  CASE ("M05")
612  s6 = 1.000_dp
613  sr6 = 1.373_dp
614  s8 = 0.595_dp
615  CASE ("M062X")
616  s6 = 1.000_dp
617  sr6 = 1.619_dp
618  s8 = 0.000_dp
619  CASE ("M06HF")
620  s6 = 1.000_dp
621  sr6 = 1.446_dp
622  s8 = 0.000_dp
623  CASE ("M06L")
624  s6 = 1.000_dp
625  sr6 = 1.581_dp
626  s8 = 0.000_dp
627  CASE ("M06N")
628  s6 = 1.000_dp
629  sr6 = 1.325_dp
630  s8 = 0.000_dp
631  CASE ("HCTH120")
632  s6 = 1.000_dp
633  sr6 = 1.221_dp
634  s8 = 1.206_dp
635  CASE ("HCTH407")
636  s6 = 1.000_dp
637  sr6 = 4.0426_dp
638  s8 = 2.7694_dp
639  CASE ("MPW2PLYP")
640  s6 = 1.000_dp
641  sr6 = 1.5527_dp
642  s8 = 0.7529_dp
643  CASE ("PKZB")
644  s6 = 1.000_dp
645  sr6 = 0.6327_dp
646  s8 = 0.000_dp
647  CASE ("PTPSS")
648  s6 = 0.750_dp
649  sr6 = 1.541_dp
650  s8 = 0.879_dp
651  CASE ("PWPB95")
652  s6 = 0.820_dp
653  sr6 = 1.557_dp
654  s8 = 0.705_dp
655  CASE ("OLYP")
656  s6 = 1.000_dp
657  sr6 = 0.806_dp
658  s8 = 1.764_dp
659  CASE ("OPBE")
660  s6 = 1.000_dp
661  sr6 = 0.837_dp
662  s8 = 2.055_dp
663  CASE ("OTPSS")
664  s6 = 1.000_dp
665  sr6 = 1.128_dp
666  s8 = 1.494_dp
667  CASE ("PBE1KCIS")
668  s6 = 1.000_dp
669  sr6 = 3.6355_dp
670  s8 = 1.7934_dp
671  CASE ("PBE38")
672  s6 = 1.000_dp
673  sr6 = 1.333_dp
674  s8 = 0.998_dp
675  CASE ("PBEH1PBE")
676  s6 = 1.000_dp
677  sr6 = 1.3719_dp
678  s8 = 1.0430_dp
679  CASE ("PBESOL")
680  s6 = 1.000_dp
681  sr6 = 1.345_dp
682  s8 = 0.612_dp
683  CASE ("REVSSB")
684  s6 = 1.000_dp
685  sr6 = 1.221_dp
686  s8 = 0.560_dp
687  CASE ("REVTPSS")
688  s6 = 1.000_dp
689  sr6 = 1.3491_dp
690  s8 = 1.3666_dp
691  CASE ("SSB")
692  s6 = 1.000_dp
693  sr6 = 1.215_dp
694  s8 = 0.663_dp
695  CASE ("B97-1")
696  s6 = 1.000_dp
697  sr6 = 3.7924_dp
698  s8 = 1.6418_dp
699  CASE ("B97-2")
700  s6 = 1.000_dp
701  sr6 = 1.7066_dp
702  s8 = 2.4661_dp
703  CASE ("B98")
704  s6 = 1.000_dp
705  sr6 = 2.6895_dp
706  s8 = 1.9078_dp
707  CASE ("BOP")
708  s6 = 1.000_dp
709  sr6 = 0.929_dp
710  s8 = 1.975_dp
711  CASE ("HISS")
712  s6 = 1.000_dp
713  sr6 = 1.3338_dp
714  s8 = 0.7615_dp
715  CASE ("HSE03")
716  s6 = 1.000_dp
717  sr6 = 1.3944_dp
718  s8 = 1.0156_dp
719  CASE ("HSE06")
720  s6 = 1.000_dp
721  sr6 = 1.129_dp
722  s8 = 0.109_dp
723  CASE ("M08HX")
724  s6 = 1.000_dp
725  sr6 = 1.6247_dp
726  s8 = 0.000_dp
727  CASE ("MN15L")
728  s6 = 1.000_dp
729  sr6 = 3.3388_dp
730  s8 = 0.000_dp
731  CASE ("MPWPW91")
732  s6 = 1.0000_dp
733  sr6 = 1.3725_dp
734  s8 = 1.9467_dp
735  CASE ("MPW1B95")
736  s6 = 1.000_dp
737  sr6 = 1.605_dp
738  s8 = 1.118_dp
739  CASE ("MPW1KCIS")
740  s6 = 1.000_dp
741  sr6 = 1.7231_dp
742  s8 = 2.2917_dp
743  CASE ("MPW1LYP")
744  s6 = 1.000_dp
745  sr6 = 2.0512_dp
746  s8 = 1.9529_dp
747  CASE ("MPW1PW91")
748  s6 = 1.000_dp
749  sr6 = 1.2892_dp
750  s8 = 1.4758_dp
751  CASE ("MPWB1K")
752  s6 = 1.000_dp
753  sr6 = 1.671_dp
754  s8 = 1.061_dp
755  CASE ("MPWKCIS1K")
756  s6 = 1.000_dp
757  sr6 = 1.4853_dp
758  s8 = 1.7553_dp
759  CASE ("O3LYP")
760  s6 = 1.000_dp
761  sr6 = 1.4060_dp
762  s8 = 1.8058_dp
763  CASE ("PW1PW")
764  s6 = 1.000_dp
765  sr6 = 1.4968_dp
766  s8 = 1.1786_dp
767  CASE ("PW91P86")
768  s6 = 1.0000_dp
769  sr6 = 2.1040_dp
770  s8 = 0.8747_dp
771  CASE ("REVPBE0")
772  s6 = 1.000_dp
773  sr6 = 0.949_dp
774  s8 = 0.792_dp
775  CASE ("REVPBE38")
776  s6 = 1.000_dp
777  sr6 = 1.021_dp
778  s8 = 0.862_dp
779  CASE ("REVTPSSh")
780  s6 = 1.000_dp
781  sr6 = 1.3224_dp
782  s8 = 1.2504_dp
783  CASE ("REVTPSS0")
784  s6 = 1.000_dp
785  sr6 = 1.2881_dp
786  s8 = 1.0649_dp
787  CASE ("TPSS1KCIS")
788  s6 = 1.000_dp
789  sr6 = 1.7729_dp
790  s8 = 2.0902_dp
791  CASE ("THCTHHYB")
792  s6 = 1.000_dp
793  sr6 = 1.5001_dp
794  s8 = 1.6302_dp
795  CASE ("RPW86PBE")
796  s6 = 1.000_dp
797  sr6 = 1.224_dp
798  s8 = 0.901_dp
799  CASE ("SCAN")
800  s6 = 1.000_dp
801  sr6 = 1.324_dp
802  s8 = 0.000_dp
803  CASE ("THCTH")
804  s6 = 1.000_dp
805  sr6 = 0.932_dp
806  s8 = 0.5662_dp
807  CASE ("XLYP")
808  s6 = 1.0000_dp
809  sr6 = 0.9384_dp
810  s8 = 0.7447_dp
811  CASE ("X3LYP")
812  s6 = 1.000_dp
813  sr6 = 1.0000_dp
814  s8 = 0.2990_dp
815  END SELECT
816 
817  END SUBROUTINE qs_scaling_dftd3
818 
819 ! **************************************************************************************************
820 !> \brief ...
821 !> \param s6 ...
822 !> \param a1 ...
823 !> \param s8 ...
824 !> \param a2 ...
825 !> \param vdw_section ...
826 ! **************************************************************************************************
827  SUBROUTINE qs_scaling_dftd3bj(s6, a1, s8, a2, vdw_section)
828  REAL(kind=dp), INTENT(inout) :: s6, a1, s8, a2
829  TYPE(section_vals_type), POINTER :: vdw_section
830 
831  CHARACTER(LEN=default_string_length) :: functional
832 
833  CALL section_vals_val_get(vdw_section, "PAIR_POTENTIAL%REFERENCE_FUNCTIONAL", c_val=functional)
834 
835  ! values for different functionals from:
836  ! http://www.thch.uni-bonn.de/tc/downloads/DFT-D3/functionalsbj.html
837  ! L. Goerigk et al. PCCP 2017, 32147-32744, SI
838  SELECT CASE (trim(functional))
839  CASE DEFAULT
840  ! unknown functional
841  cpabort("No DFT-D3(BJ) values available for this functional:"//trim(functional))
842  CASE ("B1B95")
843  s6 = 1.0000_dp
844  a1 = 0.2092_dp
845  s8 = 1.4507_dp
846  a2 = 5.5545_dp
847  CASE ("B2GPPLYP")
848  s6 = 0.5600_dp
849  a1 = 0.0000_dp
850  s8 = 0.2597_dp
851  a2 = 6.3332_dp
852  CASE ("B3PW91")
853  s6 = 1.0000_dp
854  a1 = 0.4312_dp
855  s8 = 2.8524_dp
856  a2 = 4.4693_dp
857  CASE ("BHLYP")
858  s6 = 1.0000_dp
859  a1 = 0.2793_dp
860  s8 = 1.0354_dp
861  a2 = 4.9615_dp
862  CASE ("BMK")
863  s6 = 1.0000_dp
864  a1 = 0.1940_dp
865  s8 = 2.0860_dp
866  a2 = 5.9197_dp
867  CASE ("BOP")
868  s6 = 1.0000_dp
869  a1 = 0.4870_dp
870  s8 = 3.2950_dp
871  a2 = 3.5043_dp
872  CASE ("BPBE")
873  s6 = 1.0000_dp
874  a1 = 0.4567_dp
875  s8 = 4.0728_dp
876  a2 = 4.3908_dp
877  CASE ("B97-3c")
878  s6 = 1.0000_dp
879  a1 = 0.3700_dp
880  s8 = 1.5000_dp
881  a2 = 4.1000_dp
882  CASE ("CAMB3LYP")
883  s6 = 1.0000_dp
884  a1 = 0.3708_dp
885  s8 = 2.0674_dp
886  a2 = 5.4743_dp
887  CASE ("DSDBLYP")
888  s6 = 0.5000_dp
889  a1 = 0.0000_dp
890  s8 = 0.2130_dp
891  a2 = 6.0519_dp
892  CASE ("DSDPBEP86")
893  s6 = 0.4180_dp
894  a1 = 0.0000_dp
895  s8 = 0.0000_dp
896  a2 = 5.6500_dp
897  CASE ("DSDPBEB95")
898  s6 = 0.6100_dp
899  a1 = 0.0000_dp
900  s8 = 0.0000_dp
901  a2 = 6.2000_dp
902  CASE ("LCWPBE")
903  s6 = 1.0000_dp
904  a1 = 0.3919_dp
905  s8 = 1.8541_dp
906  a2 = 5.0897_dp
907  CASE ("LCWhPBE")
908  s6 = 1.0000_dp
909  a1 = 0.2746_dp
910  s8 = 1.1908_dp
911  a2 = 5.3157_dp
912  CASE ("MPW1B95")
913  s6 = 1.0000_dp
914  a1 = 0.1955_dp
915  s8 = 1.0508_dp
916  a2 = 6.4177_dp
917  CASE ("MPW2PLYP")
918  s6 = 0.6600_dp
919  a1 = 0.4105_dp
920  s8 = 0.6223_dp
921  a2 = 5.0136_dp
922  CASE ("MPWB1K")
923  s6 = 1.0000_dp
924  a1 = 0.1474_dp
925  s8 = 0.9499_dp
926  a2 = 6.6223_dp
927  CASE ("MPWLYP")
928  s6 = 1.0000_dp
929  a1 = 0.4831_dp
930  s8 = 2.0077_dp
931  a2 = 4.5323_dp
932  CASE ("OLYP")
933  s6 = 1.0000_dp
934  a1 = 0.5299_dp
935  s8 = 2.6205_dp
936  a2 = 2.8065_dp
937  CASE ("OPBE")
938  s6 = 1.0000_dp
939  a1 = 0.5512_dp
940  s8 = 3.3816_dp
941  a2 = 2.9444_dp
942  CASE ("OTPSS")
943  s6 = 1.0000_dp
944  a1 = 0.4634_dp
945  s8 = 2.7495_dp
946  a2 = 4.3153_dp
947  CASE ("PBE38")
948  s6 = 1.0000_dp
949  a1 = 0.3995_dp
950  s8 = 1.4623_dp
951  a2 = 5.1405_dp
952  CASE ("PBEsol")
953  s6 = 1.0000_dp
954  a1 = 0.4466_dp
955  s8 = 2.9491_dp
956  a2 = 6.1742_dp
957  CASE ("PTPSS")
958  s6 = 0.7500_dp
959  a1 = 0.0000_dp
960  s8 = 0.2804_dp
961  a2 = 6.5745_dp
962  CASE ("PWB6K")
963  s6 = 1.0000_dp
964  a1 = 0.1805_dp
965  s8 = 0.9383_dp
966  a2 = 7.7627_dp
967  CASE ("revSSB")
968  s6 = 1.0000_dp
969  a1 = 0.4720_dp
970  s8 = 0.4389_dp
971  a2 = 4.0986_dp
972  CASE ("SSB")
973  s6 = 1.0000_dp
974  a1 = -0.0952_dp
975  s8 = -0.1744_dp
976  a2 = 5.2170_dp
977  CASE ("TPSSh")
978  s6 = 1.0000_dp
979  a1 = 0.4529_dp
980  s8 = 2.2382_dp
981  a2 = 4.6550_dp
982  CASE ("HCTH120")
983  s6 = 1.0000_dp
984  a1 = 0.3563_dp
985  s8 = 1.0821_dp
986  a2 = 4.3359_dp
987  CASE ("B2PLYP")
988  s6 = 0.6400_dp
989  a1 = 0.3065_dp
990  s8 = 0.9147_dp
991  a2 = 5.0570_dp
992  CASE ("B1LYP")
993  s6 = 1.0000_dp
994  a1 = 0.1986_dp
995  s8 = 2.1167_dp
996  a2 = 5.3875_dp
997  CASE ("B1P86")
998  s6 = 1.0000_dp
999  a1 = 0.4724_dp
1000  s8 = 3.5681_dp
1001  a2 = 4.9858_dp
1002  CASE ("B3LYP")
1003  s6 = 1.0000_dp
1004  a1 = 0.3981_dp
1005  s8 = 1.9889_dp
1006  a2 = 4.4211_dp
1007  CASE ("B3P86")
1008  s6 = 1.0000_dp
1009  a1 = 0.4601_dp
1010  s8 = 3.3211_dp
1011  a2 = 4.9294_dp
1012  CASE ("B97-1")
1013  s6 = 1.0000_dp
1014  a1 = 0.0000_dp
1015  s8 = 0.4814_dp
1016  a2 = 6.2279_dp
1017  CASE ("B97-2")
1018  s6 = 1.0000_dp
1019  a1 = 0.0000_dp
1020  s8 = 0.9448_dp
1021  a2 = 5.4603_dp
1022  CASE ("B97-D")
1023  s6 = 1.0000_dp
1024  a1 = 0.5545_dp
1025  s8 = 2.2609_dp
1026  a2 = 3.2297_dp
1027  CASE ("B98")
1028  s6 = 1.0000_dp
1029  a1 = 0.0000_dp
1030  s8 = 0.7086_dp
1031  a2 = 6.0672_dp
1032  CASE ("BLYP")
1033  s6 = 1.0000_dp
1034  a1 = 0.4298_dp
1035  s8 = 2.6996_dp
1036  a2 = 4.2359_dp
1037  CASE ("BP86")
1038  s6 = 1.0000_dp
1039  a1 = 0.3946_dp
1040  s8 = 3.2822_dp
1041  a2 = 4.8516_dp
1042  CASE ("DSD-BLYP")
1043  s6 = 0.5000_dp
1044  a1 = 0.0000_dp
1045  s8 = 0.2130_dp
1046  a2 = 6.0519_dp
1047  CASE ("HCTH407")
1048  s6 = 1.0000_dp
1049  a1 = 0.0000_dp
1050  s8 = 0.6490_dp
1051  a2 = 4.8162_dp
1052  CASE ("HISS")
1053  s6 = 1.0000_dp
1054  a1 = 0.0000_dp
1055  s8 = 1.6112_dp
1056  a2 = 7.3539_dp
1057  CASE ("HSE03")
1058  s6 = 1.0000_dp
1059  a1 = 0.0000_dp
1060  s8 = 1.1243_dp
1061  a2 = 6.8889_dp
1062  CASE ("HSE06")
1063  s6 = 1.0000_dp
1064  a1 = 0.3830_dp
1065  s8 = 2.3100_dp
1066  a2 = 5.6850_dp
1067  CASE ("M11")
1068  s6 = 1.0000_dp
1069  a1 = 0.0000_dp
1070  s8 = 2.8112_dp
1071  a2 = 10.1389_dp
1072  CASE ("MN12SX")
1073  s6 = 1.0000_dp
1074  a1 = 0.0983_dp
1075  s8 = 1.1674_dp
1076  a2 = 8.0259_dp
1077  CASE ("MN15")
1078  s6 = 1.0000_dp
1079  a1 = 2.0971_dp
1080  s8 = 0.7862_dp
1081  a2 = 7.5923_dp
1082  CASE ("mPWPW91")
1083  s6 = 1.0000_dp
1084  a1 = 0.3168_dp
1085  s8 = 1.7974_dp
1086  a2 = 4.7732_dp
1087  CASE ("MPW1PW91")
1088  s6 = 1.0000_dp
1089  a1 = 0.3342_dp
1090  s8 = 1.8744_dp
1091  a2 = 4.9819_dp
1092  CASE ("MPW1KCIS")
1093  s6 = 1.0000_dp
1094  a1 = 0.0576_dp
1095  s8 = 1.0893_dp
1096  a2 = 5.5314_dp
1097  CASE ("MPWKCIS1K")
1098  s6 = 1.0000_dp
1099  a1 = 0.0855_dp
1100  s8 = 1.2875_dp
1101  a2 = 5.8961_dp
1102  CASE ("N12SX")
1103  s6 = 1.0000_dp
1104  a1 = 0.3283_dp
1105  s8 = 2.4900_dp
1106  a2 = 5.7898_dp
1107  CASE ("O3LYP")
1108  s6 = 1.0000_dp
1109  a1 = 0.0963_dp
1110  s8 = 1.8171_dp
1111  a2 = 5.9940_dp
1112  CASE ("PBE0")
1113  s6 = 1.0000_dp
1114  a1 = 0.4145_dp
1115  s8 = 1.2177_dp
1116  a2 = 4.8593_dp
1117  CASE ("PBE")
1118  s6 = 1.0000_dp
1119  a1 = 0.4289_dp
1120  s8 = 0.7875_dp
1121  a2 = 4.4407_dp
1122  CASE ("PBEhPBE")
1123  s6 = 1.0000_dp
1124  a1 = 0.0000_dp
1125  s8 = 1.1152_dp
1126  a2 = 6.7184_dp
1127  CASE ("PBEh1PBE")
1128  s6 = 1.0000_dp
1129  a1 = 0.0000_dp
1130  s8 = 1.4877_dp
1131  a2 = 7.0385_dp
1132  CASE ("PBE1KCIS")
1133  s6 = 1.0000_dp
1134  a1 = 0.0000_dp
1135  s8 = 0.7688_dp
1136  a2 = 6.2794_dp
1137  CASE ("PW6B95")
1138  s6 = 1.0000_dp
1139  a1 = 0.2076_dp
1140  s8 = 0.7257_dp
1141  a2 = 6.3750_dp
1142  CASE ("PWPB95")
1143  s6 = 0.8200_dp
1144  a1 = 0.0000_dp
1145  s8 = 0.2904_dp
1146  a2 = 7.3141_dp
1147  CASE ("revPBE0")
1148  s6 = 1.0000_dp
1149  a1 = 0.4679_dp
1150  s8 = 1.7588_dp
1151  a2 = 3.7619_dp
1152  CASE ("revPBE38")
1153  s6 = 1.0000_dp
1154  a1 = 0.4309_dp
1155  s8 = 1.4760_dp
1156  a2 = 3.9446_dp
1157  CASE ("revPBE")
1158  s6 = 1.0000_dp
1159  a1 = 0.5238_dp
1160  s8 = 2.3550_dp
1161  a2 = 3.5016_dp
1162  CASE ("revTPSS")
1163  s6 = 1.0000_dp
1164  a1 = 0.4426_dp
1165  s8 = 1.4023_dp
1166  a2 = 4.4723_dp
1167  CASE ("revTPSS0")
1168  s6 = 1.0000_dp
1169  a1 = 0.2218_dp
1170  s8 = 1.6151_dp
1171  a2 = 5.7985_dp
1172  CASE ("revTPSSh")
1173  s6 = 1.0000_dp
1174  a1 = 0.2660_dp
1175  s8 = 1.4076_dp
1176  a2 = 5.3761_dp
1177  CASE ("RPBE")
1178  s6 = 1.0000_dp
1179  a1 = 0.1820_dp
1180  s8 = 0.8318_dp
1181  a2 = 4.0094_dp
1182  CASE ("RPW86PBE")
1183  s6 = 1.0000_dp
1184  a1 = 0.4613_dp
1185  s8 = 1.3845_dp
1186  a2 = 4.5062_dp
1187  CASE ("SCAN")
1188  s6 = 1.0000_dp
1189  a1 = 0.538_dp
1190  s8 = 0.0000_dp
1191  a2 = 5.420_dp
1192  CASE ("SOGGA11X")
1193  s6 = 1.0000_dp
1194  a1 = 0.1330_dp
1195  s8 = 1.1426_dp
1196  a2 = 5.7381_dp
1197  CASE ("TPSS0")
1198  s6 = 1.0000_dp
1199  a1 = 0.3768_dp
1200  s8 = 1.2576_dp
1201  a2 = 4.5865_dp
1202  CASE ("TPSS1KCIS")
1203  s6 = 1.0000_dp
1204  a1 = 0.0000_dp
1205  s8 = 1.0542_dp
1206  a2 = 6.0201_dp
1207  CASE ("TPSS")
1208  s6 = 1.0000_dp
1209  a1 = 0.4535_dp
1210  s8 = 1.9435_dp
1211  a2 = 4.4752_dp
1212  CASE ("tHCTH")
1213  s6 = 1.0000_dp
1214  a1 = 0.0000_dp
1215  s8 = 1.2626_dp
1216  a2 = 5.6162_dp
1217  CASE ("tHCTHhyb")
1218  s6 = 1.0000_dp
1219  a1 = 0.0000_dp
1220  s8 = 0.9585_dp
1221  a2 = 6.2303_dp
1222  CASE ("XLYP")
1223  s6 = 1.0000_dp
1224  a1 = 0.0809_dp
1225  s8 = 1.5669_dp
1226  a2 = 5.3166_dp
1227  CASE ("X3LYP")
1228  s6 = 1.0000_dp
1229  a1 = 0.2022_dp
1230  s8 = 1.5744_dp
1231  a2 = 5.4184_dp
1232  END SELECT
1233 
1234  END SUBROUTINE qs_scaling_dftd3bj
1235 
1236 ! **************************************************************************************************
1237 !> \brief ...
1238 !> \param qs_env ...
1239 !> \param dispersion_env ...
1240 !> \param energy ...
1241 !> \param calculate_forces ...
1242 !> \param atevdw ...
1243 ! **************************************************************************************************
1244  SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
1245 
1246  TYPE(qs_environment_type), POINTER :: qs_env
1247  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1248  REAL(kind=dp), INTENT(OUT) :: energy
1249  LOGICAL, INTENT(IN) :: calculate_forces
1250  REAL(kind=dp), DIMENSION(:), OPTIONAL :: atevdw
1251 
1252  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_dispersion_pairpot'
1253 
1254  INTEGER :: atom_a, atom_b, atom_c, atom_d, handle, hashb, hashc, i, ia, iat, iatom, icx, &
1255  icy, icz, idmp, ikind, ilist, imol, jatom, jkind, katom, kkind, kstart, latom, lkind, &
1256  max_elem, maxc, mepos, na, natom, nb, nc, nkind, num_pe, unit_nr, za, zb, zc
1257  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
1258  INTEGER, DIMENSION(3) :: cell_b, cell_c, ncell, periodic
1259  INTEGER, DIMENSION(:), POINTER :: atom_list
1260  LOGICAL :: atenergy, atex, atstress, debug, debugall, debugx, domol, exclude_pair, &
1261  floating_a, floating_b, floating_c, ghost_a, ghost_b, ghost_c, is000, use_virial
1262  LOGICAL, ALLOCATABLE, DIMENSION(:) :: dodisp, floating, ghost
1263  REAL(kind=dp) :: a1, a2, alp6, alp8, ang, c6, c8, c9, cc6ab, cc6bc, cc6ca, cnum, dc6a, dc6b, &
1264  dc8a, dc8b, dcc6aba, dcc6abb, dcc6bcb, dcc6bcc, dcc6caa, dcc6cac, dd, de6, de8, de91, &
1265  de921, de922, dea, devdw, dfdab6, dfdab8, dfdabc, dfdmp, dr, drk, e6, e6tot, e8, e8tot, &
1266  e9, e9tot, elrc6, elrc8, elrc9, eps_cn, er, esrb, evdw, f0ab, fac, fac0, fdab6, fdab8, &
1267  fdabc, fdmp, gnorm, gsrb, kgc8, nab, nabc, r0, r0ab, r2ab, r2abc, r2bc, r2ca, r6, r8, &
1268  rabc, rc2, rcc, rcut, rs6, rs8, s1, s2, s3, s6, s8, s8i, s9, srbe, ssrb, t1srb, t2srb, xp
1269  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: atom2mol, c6d2, cnkind, cnumbers, &
1270  cnumfix, radd2
1271  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rcpbc
1272  REAL(kind=dp), DIMENSION(3) :: fdij, fdik, ra, rab, rb, rb0, rbc, rc, &
1273  rc0, rca, rij, rik, sab_max
1274  REAL(kind=dp), DIMENSION(3, 3) :: dvirial, pv_loc, pv_virial_thread
1275  REAL(kind=dp), DIMENSION(:), POINTER :: atener
1276  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: atstr
1277  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1278  TYPE(atprop_type), POINTER :: atprop
1279  TYPE(cell_type), POINTER :: cell
1280  TYPE(cp_logger_type), POINTER :: logger
1281  TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
1282  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1283  TYPE(mp_para_env_type), POINTER :: para_env
1284  TYPE(neighbor_list_iterator_p_type), &
1285  DIMENSION(:), POINTER :: nl_iterator
1286  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1287  POINTER :: sab_vdw
1288  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1289  TYPE(qs_atom_dispersion_type), POINTER :: disp_a, disp_b, disp_c
1290  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1291  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1292  TYPE(virial_type), POINTER :: virial
1293 
1294  energy = 0._dp
1295  ! make valgrind happy
1296  use_virial = .false.
1297 
1298  IF (dispersion_env%type .NE. xc_vdw_fun_pairpot) THEN
1299  RETURN
1300  END IF
1301 
1302  CALL timeset(routinen, handle)
1303 
1304  NULLIFY (atomic_kind_set, qs_kind_set, sab_vdw)
1305 
1306  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1307  cell=cell, virial=virial, para_env=para_env, atprop=atprop)
1308 
1309  debugx = dispersion_env%verbose
1310  debugall = debugx
1311  debug = debugx .AND. para_env%is_source()
1312  nkind = SIZE(atomic_kind_set)
1313 
1314  NULLIFY (logger)
1315  logger => cp_get_default_logger()
1316  IF (ASSOCIATED(dispersion_env%dftd_section)) THEN
1317  unit_nr = cp_print_key_unit_nr(logger, dispersion_env%dftd_section, "PRINT_DFTD", &
1318  extension=".dftd")
1319  ELSE
1320  unit_nr = -1
1321  END IF
1322 
1323  ! atomic energy and stress arrays
1324  atenergy = atprop%energy
1325  IF (atenergy) THEN
1326  NULLIFY (particle_set)
1327  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
1328  natom = SIZE(particle_set)
1329  CALL atprop_array_init(atprop%atevdw, natom)
1330  atener => atprop%atevdw
1331  END IF
1332  atstress = atprop%stress
1333  atstr => atprop%atstress
1334  ! external atomic energy
1335  atex = .false.
1336  IF (PRESENT(atevdw)) THEN
1337  atex = .true.
1338  END IF
1339 
1340  IF (unit_nr > 0) THEN
1341  WRITE (unit_nr, *)
1342  WRITE (unit_nr, *) " Pair potential vdW calculation"
1343  IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
1344  WRITE (unit_nr, *) " Dispersion potential type: DFT-D2"
1345  ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
1346  WRITE (unit_nr, *) " Dispersion potential type: DFT-D3"
1347  ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
1348  WRITE (unit_nr, *) " Dispersion potential type: DFT-D3(BJ)"
1349  END IF
1350  END IF
1351 
1352  NULLIFY (particle_set)
1353  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
1354  natom = SIZE(particle_set)
1355  IF (calculate_forces .OR. debugall) THEN
1356  NULLIFY (force)
1357  CALL get_qs_env(qs_env=qs_env, force=force)
1358  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
1359  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1360  IF (use_virial .AND. debugall) THEN
1361  dvirial = virial%pv_virial
1362  END IF
1363  IF (use_virial) THEN
1364  pv_loc = virial%pv_virial
1365  END IF
1366  ELSE IF ((dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
1367  dispersion_env%pp_type == vdw_pairpot_dftd3bj) .AND. dispersion_env%doabc) THEN
1368  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
1369  END IF
1370 
1371  ALLOCATE (dodisp(nkind), ghost(nkind), floating(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
1372  DO ikind = 1, nkind
1373  CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
1374  CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
1375  dodisp(ikind) = disp_a%defined
1376  ghost(ikind) = ghost_a
1377  floating(ikind) = floating_a
1378  atomnumber(ikind) = za
1379  c6d2(ikind) = disp_a%c6
1380  radd2(ikind) = disp_a%vdw_radii
1381  END DO
1382 
1383  ALLOCATE (rcpbc(3, natom))
1384  DO iatom = 1, natom
1385  rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
1386  END DO
1387 
1388  rcut = 2._dp*dispersion_env%rc_disp
1389  rc2 = rcut*rcut
1390  IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
1391  s6 = dispersion_env%scaling
1392  dd = dispersion_env%exp_pre
1393  IF (unit_nr > 0) THEN
1394  WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
1395  WRITE (unit_nr, *) " Exponential prefactor ", dd
1396  END IF
1397  domol = .false.
1398  ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
1399  dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
1400  maxc = SIZE(dispersion_env%c6ab, 3)
1401  max_elem = SIZE(dispersion_env%c6ab, 1)
1402  alp6 = dispersion_env%alp
1403  alp8 = alp6 + 2._dp
1404  s6 = dispersion_env%s6
1405  s8 = dispersion_env%s8
1406  s9 = dispersion_env%s6
1407  rs6 = dispersion_env%sr6
1408  rs8 = 1._dp
1409  a1 = dispersion_env%a1
1410  a2 = dispersion_env%a2
1411  eps_cn = dispersion_env%eps_cn
1412  e6tot = 0._dp
1413  e8tot = 0._dp
1414  e9tot = 0._dp
1415  esrb = 0._dp
1416  domol = dispersion_env%domol
1417  ! molecule correction
1418  kgc8 = dispersion_env%kgc8
1419  IF (domol) THEN
1420  CALL get_qs_env(qs_env=qs_env, molecule_set=molecule_set)
1421  ALLOCATE (atom2mol(natom))
1422  DO imol = 1, SIZE(molecule_set)
1423  DO iat = molecule_set(imol)%first_atom, molecule_set(imol)%last_atom
1424  atom2mol(iat) = imol
1425  END DO
1426  END DO
1427  END IF
1428  ! damping type
1429  idmp = 0
1430  IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
1431  idmp = 1
1432  ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
1433  idmp = 2
1434  END IF
1435  ! SRB parameters
1436  ssrb = dispersion_env%srb_params(1)
1437  gsrb = dispersion_env%srb_params(2)
1438  t1srb = dispersion_env%srb_params(3)
1439  t2srb = dispersion_env%srb_params(4)
1440 
1441  IF (unit_nr > 0) THEN
1442  WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
1443  WRITE (unit_nr, *) " Scaling parameter (s8) ", s8
1444  IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
1445  WRITE (unit_nr, *) " Zero Damping parameter (sr6)", rs6
1446  WRITE (unit_nr, *) " Zero Damping parameter (sr8)", rs8
1447  ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
1448  WRITE (unit_nr, *) " BJ Damping parameter (a1) ", a1
1449  WRITE (unit_nr, *) " BJ Damping parameter (a2) ", a2
1450  END IF
1451  WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
1452  IF (dispersion_env%lrc) THEN
1453  WRITE (unit_nr, *) " Apply a long range correction"
1454  END IF
1455  IF (dispersion_env%srb) THEN
1456  WRITE (unit_nr, *) " Apply a short range bond correction"
1457  WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
1458  END IF
1459  IF (domol) THEN
1460  WRITE (unit_nr, *) " Inter-molecule scaling parameter (s8) ", kgc8
1461  END IF
1462  END IF
1463  ! Calculate coordination numbers
1464  NULLIFY (particle_set)
1465  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
1466  natom = SIZE(particle_set)
1467  ALLOCATE (cnumbers(natom))
1468  cnumbers = 0._dp
1469 
1470  IF (calculate_forces .OR. debugall) THEN
1471  ALLOCATE (dcnum(natom))
1472  dcnum(:)%neighbors = 0
1473  DO iatom = 1, natom
1474  ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
1475  END DO
1476  ELSE
1477  ALLOCATE (dcnum(1))
1478  END IF
1479 
1480  CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, &
1481  calculate_forces, debugall)
1482 
1483  CALL para_env%sum(cnumbers)
1484  ! for parallel runs we have to update dcnum on all processors
1485  IF (calculate_forces .OR. debugall) THEN
1486  CALL dcnum_distribute(dcnum, para_env)
1487  IF (unit_nr > 0 .AND. SIZE(dcnum) > 0) THEN
1488  WRITE (unit_nr, *)
1489  WRITE (unit_nr, *) " ATOM Coordination Neighbors"
1490  DO i = 1, natom
1491  WRITE (unit_nr, "(I6,F20.10,I12)") i, cnumbers(i), dcnum(i)%neighbors
1492  END DO
1493  WRITE (unit_nr, *)
1494  END IF
1495  END IF
1496  END IF
1497 
1498  nab = 0._dp
1499  nabc = 0._dp
1500  IF (dispersion_env%doabc) THEN
1501  rcc = 2._dp*dispersion_env%rc_disp
1502  CALL get_cell(cell=cell, periodic=periodic)
1503  sab_max(1) = rcc/plane_distance(1, 0, 0, cell)
1504  sab_max(2) = rcc/plane_distance(0, 1, 0, cell)
1505  sab_max(3) = rcc/plane_distance(0, 0, 1, cell)
1506  ncell(:) = (int(sab_max(:)) + 1)*periodic(:)
1507  IF (unit_nr > 0) THEN
1508  WRITE (unit_nr, *) " Calculate C9 Terms"
1509  WRITE (unit_nr, "(A,T20,I3,A,I3)") " Search in cells ", -ncell(1), ":", ncell(1)
1510  WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(2), ":", ncell(2)
1511  WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(3), ":", ncell(3)
1512  WRITE (unit_nr, *)
1513  END IF
1514  IF (dispersion_env%c9cnst) THEN
1515  IF (unit_nr > 0) WRITE (unit_nr, *) " Use reference coordination numbers for C9 term"
1516  ALLOCATE (cnumfix(natom))
1517  cnumfix = 0._dp
1518  ! first use the default values
1519  DO iatom = 1, natom
1520  ikind = kind_of(iatom)
1521  CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
1522  cnumfix(iatom) = dispersion_env%cn(za)
1523  END DO
1524  ! now check for changes from default
1525  IF (ASSOCIATED(dispersion_env%cnkind)) THEN
1526  DO i = 1, SIZE(dispersion_env%cnkind)
1527  ikind = dispersion_env%cnkind(i)%kind
1528  cnum = dispersion_env%cnkind(i)%cnum
1529  cpassert(ikind <= nkind)
1530  cpassert(ikind > 0)
1531  CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, atom_list=atom_list)
1532  DO ia = 1, na
1533  iatom = atom_list(ia)
1534  cnumfix(iatom) = cnum
1535  END DO
1536  END DO
1537  END IF
1538  IF (ASSOCIATED(dispersion_env%cnlist)) THEN
1539  DO i = 1, SIZE(dispersion_env%cnlist)
1540  DO ilist = 1, dispersion_env%cnlist(i)%natom
1541  iatom = dispersion_env%cnlist(i)%atom(ilist)
1542  cnumfix(iatom) = dispersion_env%cnlist(i)%cnum
1543  END DO
1544  END DO
1545  END IF
1546  IF (unit_nr > 0) THEN
1547  DO i = 1, natom
1548  IF (abs(cnumbers(i) - cnumfix(i)) > 0.5_dp) THEN
1549  WRITE (unit_nr, "(A,T20,A,I6,T41,2F10.3)") " Difference in CN ", "Atom:", &
1550  i, cnumbers(i), cnumfix(i)
1551  END IF
1552  END DO
1553  WRITE (unit_nr, *)
1554  END IF
1555  END IF
1556  END IF
1557 
1558  evdw = 0._dp
1559  sab_vdw => dispersion_env%sab_vdw
1560  nkind = SIZE(atomic_kind_set)
1561 
1562  num_pe = 1
1563 
1564  pv_virial_thread(:, :) = 0._dp
1565 
1566  CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
1567 
1568  mepos = 0
1569  DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
1570  CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
1571 
1572  IF (ghost(ikind) .OR. ghost(jkind) .OR. floating(ikind) .OR. floating(jkind)) cycle
1573 
1574  IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) cycle
1575 
1576  za = atomnumber(ikind)
1577  zb = atomnumber(jkind)
1578  ! vdW potential
1579  dr = sqrt(sum(rij(:)**2))
1580  IF (dr <= rcut) THEN
1581  nab = nab + 1._dp
1582  fac = 1._dp
1583  IF (iatom == jatom) fac = 0.5_dp
1584  IF (disp_a%type == dftd2_pp .AND. dr > 0.001_dp) THEN
1585  c6 = sqrt(c6d2(ikind)*c6d2(jkind))
1586  rcc = radd2(ikind) + radd2(jkind)
1587  er = exp(-dd*(dr/rcc - 1._dp))
1588  fdmp = 1._dp/(1._dp + er)
1589  xp = s6*c6/dr**6
1590  evdw = evdw - xp*fdmp*fac
1591  IF (calculate_forces) THEN
1592  dfdmp = dd/rcc*er*fdmp*fdmp
1593  devdw = -xp*(-6._dp*fdmp/dr + dfdmp)
1594  fdij(:) = devdw*rij(:)/dr*fac
1595  atom_a = atom_of_kind(iatom)
1596  atom_b = atom_of_kind(jatom)
1597  force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
1598  force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
1599  IF (use_virial) THEN
1600  CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
1601  END IF
1602  IF (atstress) THEN
1603  CALL virial_pair_force(atstr(:, :, iatom), -0.5_dp, fdij, rij)
1604  CALL virial_pair_force(atstr(:, :, jatom), -0.5_dp, fdij, rij)
1605  END IF
1606  END IF
1607  IF (atenergy) THEN
1608  atener(iatom) = atener(iatom) - 0.5_dp*xp*fdmp*fac
1609  atener(jatom) = atener(jatom) - 0.5_dp*xp*fdmp*fac
1610  END IF
1611  IF (atex) THEN
1612  atevdw(iatom) = atevdw(iatom) - 0.5_dp*xp*fdmp*fac
1613  atevdw(jatom) = atevdw(jatom) - 0.5_dp*xp*fdmp*fac
1614  END IF
1615  ELSE IF (disp_a%type == dftd3_pp .AND. dr > 0.001_dp) THEN
1616  IF (dispersion_env%nd3_exclude_pair > 0) THEN
1617  CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, exclude=exclude_pair)
1618  IF (exclude_pair) cycle
1619  END IF
1620  ! C6 coefficient
1621  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
1622  cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, c6, dc6a, dc6b)
1623  c8 = 3.0d0*c6*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
1624  dc8a = 3.0d0*dc6a*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
1625  dc8b = 3.0d0*dc6b*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
1626  r6 = dr**6
1627  r8 = r6*dr*dr
1628  s8i = s8
1629  IF (domol) THEN
1630  IF (atom2mol(iatom) /= atom2mol(jatom)) THEN
1631  s8i = kgc8
1632  END IF
1633  END IF
1634  ! damping
1635  IF (idmp == 1) THEN
1636  ! zero
1637  CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs6, alp6, rcut, fdab6, dfdab6)
1638  CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs8, alp8, rcut, fdab8, dfdab8)
1639  e6 = s6*fac*c6*fdab6/r6
1640  e8 = s8i*fac*c8*fdab8/r8
1641  ELSE IF (idmp == 2) THEN
1642  ! BJ
1643  r0ab = sqrt(3.0d0*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb))
1644  f0ab = a1*r0ab + a2
1645  fdab6 = 1.0_dp/(r6 + f0ab**6)
1646  fdab8 = 1.0_dp/(r8 + f0ab**8)
1647  e6 = s6*fac*c6*fdab6
1648  e8 = s8i*fac*c8*fdab8
1649  ELSE
1650  cpabort("Unknown DFT-D3 damping function:")
1651  END IF
1652  IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
1653  srbe = ssrb*(real((za*zb), kind=dp))**t1srb*exp(-gsrb*dr*dispersion_env%r0ab(za, zb)**t2srb)
1654  esrb = esrb + srbe
1655  evdw = evdw - srbe
1656  ELSE
1657  srbe = 0.0_dp
1658  END IF
1659  evdw = evdw - e6 - e8
1660  e6tot = e6tot - e6
1661  e8tot = e8tot - e8
1662  IF (atenergy) THEN
1663  atener(iatom) = atener(iatom) - 0.5_dp*(e6 + e8 + srbe)
1664  atener(jatom) = atener(jatom) - 0.5_dp*(e6 + e8 + srbe)
1665  END IF
1666  IF (atex) THEN
1667  atevdw(iatom) = atevdw(iatom) - 0.5_dp*(e6 + e8 + srbe)
1668  atevdw(jatom) = atevdw(jatom) - 0.5_dp*(e6 + e8 + srbe)
1669  END IF
1670  IF (calculate_forces) THEN
1671  ! damping
1672  IF (idmp == 1) THEN
1673  ! zero
1674  de6 = -s6*c6/r6*(dfdab6 - 6._dp*fdab6/dr)
1675  de8 = -s8i*c8/r8*(dfdab8 - 8._dp*fdab8/dr)
1676  ELSE IF (idmp == 2) THEN
1677  ! BJ
1678  de6 = s6*c6*fdab6*fdab6*6.0_dp*dr**5
1679  de8 = s8i*c8*fdab8*fdab8*8.0_dp*dr**7
1680  ELSE
1681  cpabort("Unknown DFT-D3 damping function:")
1682  END IF
1683  fdij(:) = (de6 + de8)*rij(:)/dr*fac
1684  IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
1685  fdij(:) = fdij(:) + srbe*gsrb*dispersion_env%r0ab(za, zb)**t2srb*rij(:)/dr
1686  END IF
1687  atom_a = atom_of_kind(iatom)
1688  atom_b = atom_of_kind(jatom)
1689  force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
1690  force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
1691  IF (use_virial) THEN
1692  CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
1693  END IF
1694  IF (atstress) THEN
1695  CALL virial_pair_force(atstr(:, :, iatom), -0.5_dp, fdij, rij)
1696  CALL virial_pair_force(atstr(:, :, jatom), -0.5_dp, fdij, rij)
1697  END IF
1698  ! forces from the r-dependence of the coordination numbers
1699  IF (idmp == 1) THEN
1700  ! zero
1701  dc6a = -s6*fac*dc6a*fdab6/r6
1702  dc6b = -s6*fac*dc6b*fdab6/r6
1703  dc8a = -s8i*fac*dc8a*fdab8/r8
1704  dc8b = -s8i*fac*dc8b*fdab8/r8
1705  ELSE IF (idmp == 2) THEN
1706  ! BJ
1707  dc6a = -s6*fac*dc6a*fdab6
1708  dc6b = -s6*fac*dc6b*fdab6
1709  dc8a = -s8i*fac*dc8a*fdab8
1710  dc8b = -s8i*fac*dc8b*fdab8
1711  ELSE
1712  cpabort("Unknown DFT-D3 damping function:")
1713  END IF
1714  DO i = 1, dcnum(iatom)%neighbors
1715  katom = dcnum(iatom)%nlist(i)
1716  kkind = kind_of(katom)
1717  rik = dcnum(iatom)%rik(:, i)
1718  drk = sqrt(sum(rik(:)**2))
1719  fdik(:) = (dc6a + dc8a)*dcnum(iatom)%dvals(i)*rik(:)/drk
1720  atom_c = atom_of_kind(katom)
1721  force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
1722  force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
1723  IF (use_virial) THEN
1724  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1725  END IF
1726  IF (atstress) THEN
1727  CALL virial_pair_force(atstr(:, :, iatom), -0.5_dp, fdik, rik)
1728  CALL virial_pair_force(atstr(:, :, katom), -0.5_dp, fdik, rik)
1729  END IF
1730  END DO
1731  DO i = 1, dcnum(jatom)%neighbors
1732  katom = dcnum(jatom)%nlist(i)
1733  kkind = kind_of(katom)
1734  rik = dcnum(jatom)%rik(:, i)
1735  drk = sqrt(sum(rik(:)**2))
1736  fdik(:) = (dc6b + dc8b)*dcnum(jatom)%dvals(i)*rik(:)/drk
1737  atom_c = atom_of_kind(katom)
1738  force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
1739  force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
1740  IF (use_virial) THEN
1741  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1742  END IF
1743  IF (atstress) THEN
1744  CALL virial_pair_force(atstr(:, :, jatom), -0.5_dp, fdik, rik)
1745  CALL virial_pair_force(atstr(:, :, katom), -0.5_dp, fdik, rik)
1746  END IF
1747  END DO
1748  END IF
1749  IF (dispersion_env%doabc) THEN
1750  CALL get_iterator_info(nl_iterator, cell=cell_b)
1751  hashb = cellhash(cell_b, ncell)
1752  is000 = (all(cell_b == 0))
1753  rb0(:) = matmul(cell%hmat, cell_b)
1754  ra(:) = pbc(particle_set(iatom)%r(:), cell)
1755  rb(:) = pbc(particle_set(jatom)%r(:), cell) + rb0
1756  r2ab = sum((ra - rb)**2)
1757  DO icx = -ncell(1), ncell(1)
1758  DO icy = -ncell(2), ncell(2)
1759  DO icz = -ncell(3), ncell(3)
1760  cell_c(1) = icx
1761  cell_c(2) = icy
1762  cell_c(3) = icz
1763  hashc = cellhash(cell_c, ncell)
1764  IF (is000 .AND. (all(cell_c == 0))) THEN
1765  ! CASE 1: all atoms in (000), use only ordered triples
1766  kstart = max(jatom + 1, iatom + 1)
1767  fac0 = 1.0_dp
1768  ELSE IF (is000) THEN
1769  ! CASE 2: AB in (000), C in other cell
1770  ! This case covers also all instances with BC in same
1771  ! cell not (000)
1772  kstart = 1
1773  fac0 = 1.0_dp
1774  ELSE
1775  ! These are case 2 again, cycle
1776  IF (hashc == hashb) cycle
1777  IF (all(cell_c == 0)) cycle
1778  ! CASE 3: A in (000) and B and C in different cells
1779  kstart = 1
1780  fac0 = 1.0_dp/3.0_dp
1781  END IF
1782  rc0 = matmul(cell%hmat, cell_c)
1783  DO katom = kstart, natom
1784  kkind = kind_of(katom)
1785  IF (ghost(kkind) .OR. floating(kkind) .OR. .NOT. dodisp(kkind)) cycle
1786  rc(:) = rcpbc(:, katom) + rc0(:)
1787  r2bc = sum((rb - rc)**2)
1788  IF (r2bc >= rc2) cycle
1789  r2ca = sum((rc - ra)**2)
1790  IF (r2ca >= rc2) cycle
1791  r2abc = r2ab*r2bc*r2ca
1792  IF (r2abc <= 0.001_dp) cycle
1793  IF (dispersion_env%nd3_exclude_pair > 0) THEN
1794  CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
1795  kkind, exclude_pair)
1796  IF (exclude_pair) cycle
1797  END IF
1798  ! this is an empirical scaling
1799  IF (r2abc <= 0.01*rc2*rc2*rc2) THEN
1800  kkind = kind_of(katom)
1801  atom_c = atom_of_kind(katom)
1802  zc = atomnumber(kkind)
1803  ! avoid double counting!
1804  fac = 1._dp
1805  IF (iatom == jatom .OR. iatom == katom .OR. jatom == katom) fac = 0.5_dp
1806  IF (iatom == jatom .AND. iatom == katom) fac = 1._dp/3._dp
1807  fac = fac*fac0
1808  nabc = nabc + 1._dp
1809  IF (dispersion_env%c9cnst) THEN
1810  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
1811  cnumfix(iatom), cnumfix(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
1812  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
1813  cnumfix(jatom), cnumfix(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
1814  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
1815  cnumfix(katom), cnumfix(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
1816  ELSE
1817  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
1818  cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
1819  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
1820  cnumbers(jatom), cnumbers(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
1821  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
1822  cnumbers(katom), cnumbers(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
1823  END IF
1824  c9 = -sqrt(cc6ab*cc6bc*cc6ca)
1825  rabc = r2abc**(1._dp/6._dp)
1826  r0 = (dispersion_env%r0ab(za, zb)*dispersion_env%r0ab(zb, zc)* &
1827  dispersion_env%r0ab(zc, za))**(1._dp/3._dp)
1828  ! bug fixed 3.10.2017
1829  ! correct value from alp6=14 to 16 as used in original paper
1830  ! CALL damping_d3(rabc, r0, 4._dp/3._dp, alp6, rcut, fdabc, dfdabc)
1831  CALL damping_d3(rabc, r0, 4._dp/3._dp, 16.0_dp, rcut, fdabc, dfdabc)
1832  s1 = r2ab + r2bc - r2ca
1833  s2 = r2ab + r2ca - r2bc
1834  s3 = r2ca + r2bc - r2ab
1835  ang = 0.375_dp*s1*s2*s3/r2abc + 1.0_dp
1836 
1837  e9 = s9*fac*fdabc*c9*ang/r2abc**1.50d0
1838  evdw = evdw - e9
1839  e9tot = e9tot - e9
1840  IF (atenergy) THEN
1841  atener(iatom) = atener(iatom) - e9/3._dp
1842  atener(jatom) = atener(jatom) - e9/3._dp
1843  atener(katom) = atener(katom) - e9/3._dp
1844  END IF
1845  IF (atex) THEN
1846  atevdw(iatom) = atevdw(iatom) - e9/3._dp
1847  atevdw(jatom) = atevdw(jatom) - e9/3._dp
1848  atevdw(katom) = atevdw(katom) - e9/3._dp
1849  END IF
1850 
1851  IF (calculate_forces) THEN
1852  rab = rb - ra; rbc = rc - rb; rca = ra - rc
1853  de91 = s9*fac*dfdabc*c9*ang/r2abc**1.50d0
1854  de91 = -de91/3._dp*rabc + 3._dp*e9
1855  dea = s9*fac*fdabc*c9/r2abc**2.50d0*0.75_dp
1856  fdij(:) = de91*rab(:)/r2ab
1857  fdij(:) = fdij(:) + dea*s1*s2*s3*rab(:)/r2ab
1858  fdij(:) = fdij(:) - dea*(s2*s3 + s1*s3 - s1*s2)*rab(:)
1859  force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
1860  force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
1861  IF (use_virial) THEN
1862  CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rab)
1863  END IF
1864  IF (atstress) THEN
1865  CALL virial_pair_force(atstr(:, :, iatom), -0.5_dp, fdij, rab)
1866  CALL virial_pair_force(atstr(:, :, jatom), -0.5_dp, fdij, rab)
1867  END IF
1868  fdij(:) = de91*rbc(:)/r2bc
1869  fdij(:) = fdij(:) + dea*s1*s2*s3*rbc(:)/r2bc
1870  fdij(:) = fdij(:) - dea*(s2*s3 - s1*s3 + s1*s2)*rbc(:)
1871  force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdij(:)
1872  force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdij(:)
1873  IF (use_virial) THEN
1874  CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rbc)
1875  END IF
1876  IF (atstress) THEN
1877  CALL virial_pair_force(atstr(:, :, jatom), -0.5_dp, fdij, rbc)
1878  CALL virial_pair_force(atstr(:, :, katom), -0.5_dp, fdij, rbc)
1879  END IF
1880  fdij(:) = de91*rca(:)/r2ca
1881  fdij(:) = fdij(:) + dea*s1*s2*s3*rca(:)/r2ca
1882  fdij(:) = fdij(:) - dea*(-s2*s3 + s1*s3 + s1*s2)*rca(:)
1883  force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdij(:)
1884  force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) + fdij(:)
1885  IF (use_virial) THEN
1886  CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rca)
1887  END IF
1888  IF (atstress) THEN
1889  CALL virial_pair_force(atstr(:, :, iatom), -0.5_dp, fdij, rca)
1890  CALL virial_pair_force(atstr(:, :, katom), -0.5_dp, fdij, rca)
1891  END IF
1892 
1893  IF (.NOT. dispersion_env%c9cnst) THEN
1894  ! forces from the r-dependence of the coordination numbers
1895  ! atomic stress not implemented
1896 
1897  de91 = 0.5_dp*e9/cc6ab
1898  de921 = de91*dcc6aba
1899  de922 = de91*dcc6abb
1900  DO i = 1, dcnum(iatom)%neighbors
1901  latom = dcnum(iatom)%nlist(i)
1902  lkind = kind_of(latom)
1903  rik(1) = dcnum(iatom)%rik(1, i)
1904  rik(2) = dcnum(iatom)%rik(2, i)
1905  rik(3) = dcnum(iatom)%rik(3, i)
1906  drk = sqrt(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
1907  fdik(:) = -de921*dcnum(iatom)%dvals(i)*rik(:)/drk
1908  atom_d = atom_of_kind(latom)
1909  force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
1910  force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
1911  IF (use_virial) THEN
1912  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1913  END IF
1914  END DO
1915  DO i = 1, dcnum(jatom)%neighbors
1916  latom = dcnum(jatom)%nlist(i)
1917  lkind = kind_of(latom)
1918  rik(1) = dcnum(jatom)%rik(1, i)
1919  rik(2) = dcnum(jatom)%rik(2, i)
1920  rik(3) = dcnum(jatom)%rik(3, i)
1921  drk = sqrt(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
1922  fdik(:) = -de922*dcnum(jatom)%dvals(i)*rik(:)/drk
1923  atom_d = atom_of_kind(latom)
1924  force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
1925  force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
1926  IF (use_virial) THEN
1927  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1928  END IF
1929  END DO
1930 
1931  de91 = 0.5_dp*e9/cc6bc
1932  de921 = de91*dcc6bcb
1933  de922 = de91*dcc6bcc
1934  DO i = 1, dcnum(jatom)%neighbors
1935  latom = dcnum(jatom)%nlist(i)
1936  lkind = kind_of(latom)
1937  rik(1) = dcnum(jatom)%rik(1, i)
1938  rik(2) = dcnum(jatom)%rik(2, i)
1939  rik(3) = dcnum(jatom)%rik(3, i)
1940  drk = sqrt(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
1941  fdik(:) = -de921*dcnum(jatom)%dvals(i)*rik(:)/drk
1942  atom_d = atom_of_kind(latom)
1943  force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
1944  force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
1945  IF (use_virial) THEN
1946  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1947  END IF
1948  END DO
1949  DO i = 1, dcnum(katom)%neighbors
1950  latom = dcnum(katom)%nlist(i)
1951  lkind = kind_of(latom)
1952  rik(1) = dcnum(katom)%rik(1, i)
1953  rik(2) = dcnum(katom)%rik(2, i)
1954  rik(3) = dcnum(katom)%rik(3, i)
1955  drk = sqrt(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
1956  fdik(:) = -de922*dcnum(katom)%dvals(i)*rik(:)/drk
1957  atom_d = atom_of_kind(latom)
1958  force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
1959  force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
1960  IF (use_virial) THEN
1961  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1962  END IF
1963  END DO
1964 
1965  de91 = 0.5_dp*e9/cc6ca
1966  de921 = de91*dcc6cac
1967  de922 = de91*dcc6caa
1968  DO i = 1, dcnum(katom)%neighbors
1969  latom = dcnum(katom)%nlist(i)
1970  lkind = kind_of(latom)
1971  rik(1) = dcnum(katom)%rik(1, i)
1972  rik(2) = dcnum(katom)%rik(2, i)
1973  rik(3) = dcnum(katom)%rik(3, i)
1974  drk = sqrt(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
1975  fdik(:) = -de921*dcnum(katom)%dvals(i)*rik(:)/drk
1976  atom_d = atom_of_kind(latom)
1977  force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
1978  force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
1979  IF (use_virial) THEN
1980  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1981  END IF
1982  END DO
1983  DO i = 1, dcnum(iatom)%neighbors
1984  latom = dcnum(iatom)%nlist(i)
1985  lkind = kind_of(latom)
1986  rik(1) = dcnum(iatom)%rik(1, i)
1987  rik(2) = dcnum(iatom)%rik(2, i)
1988  rik(3) = dcnum(iatom)%rik(3, i)
1989  drk = sqrt(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
1990  fdik(:) = -de922*dcnum(iatom)%dvals(i)*rik(:)/drk
1991  atom_d = atom_of_kind(latom)
1992  force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
1993  force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
1994  IF (use_virial) THEN
1995  CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
1996  END IF
1997  END DO
1998  END IF
1999 
2000  END IF
2001 
2002  END IF
2003  END DO
2004  END DO
2005  END DO
2006  END DO
2007 
2008  END IF
2009  END IF
2010  END IF
2011  END DO
2012 
2013  virial%pv_virial = virial%pv_virial + pv_virial_thread
2014 
2015  CALL neighbor_list_iterator_release(nl_iterator)
2016 
2017  elrc6 = 0._dp
2018  elrc8 = 0._dp
2019  elrc9 = 0._dp
2020  IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
2021  dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
2022  ! Long range correction (atomic contributions not implemented)
2023  IF (dispersion_env%lrc) THEN
2024  ALLOCATE (cnkind(nkind))
2025  cnkind = 0._dp
2026  ! first use the default values
2027  DO ikind = 1, nkind
2028  CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
2029  cnkind(ikind) = dispersion_env%cn(za)
2030  END DO
2031  ! now check for changes from default
2032  IF (ASSOCIATED(dispersion_env%cnkind)) THEN
2033  DO i = 1, SIZE(dispersion_env%cnkind)
2034  ikind = dispersion_env%cnkind(i)%kind
2035  cnkind(ikind) = dispersion_env%cnkind(i)%cnum
2036  END DO
2037  END IF
2038  DO ikind = 1, nkind
2039  CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, z=za)
2040  CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
2041  IF (.NOT. disp_a%defined .OR. ghost_a .OR. floating_a) cycle
2042  DO jkind = 1, nkind
2043  CALL get_atomic_kind(atomic_kind_set(jkind), natom=nb, z=zb)
2044  CALL get_qs_kind(qs_kind_set(jkind), dispersion=disp_b, ghost=ghost_b, floating=floating_b)
2045  IF (.NOT. disp_b%defined .OR. ghost_b .OR. floating_b) cycle
2046  IF (dispersion_env%nd3_exclude_pair > 0) THEN
2047  CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, exclude=exclude_pair)
2048  IF (exclude_pair) cycle
2049  END IF
2050  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
2051  cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
2052  elrc6 = elrc6 - s6*twopi*real(na*nb, kind=dp)*cc6ab/(3._dp*rcut**3*cell%deth)
2053  c8 = 3.0d0*cc6ab*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
2054  elrc8 = elrc8 - s8*twopi*real(na*nb, kind=dp)*c8/(5._dp*rcut**5*cell%deth)
2055  IF (dispersion_env%doabc) THEN
2056  DO kkind = 1, nkind
2057  CALL get_atomic_kind(atomic_kind_set(kkind), natom=nc, z=zc)
2058  CALL get_qs_kind(qs_kind_set(kkind), dispersion=disp_c, ghost=ghost_c, floating=floating_c)
2059  IF (.NOT. disp_c%defined .OR. ghost_c .OR. floating_c) cycle
2060  IF (dispersion_env%nd3_exclude_pair > 0) THEN
2061  CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, kkind, exclude_pair)
2062  IF (exclude_pair) cycle
2063  END IF
2064  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
2065  cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
2066  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
2067  cnkind(kkind), cnkind(ikind), dispersion_env%k3, cc6ca, dcc6aba, dcc6abb)
2068  CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
2069  cnkind(jkind), cnkind(kkind), dispersion_env%k3, cc6bc, dcc6aba, dcc6abb)
2070  c9 = -sqrt(cc6ab*cc6bc*cc6ca)
2071  elrc9 = elrc9 - s9*64._dp*twopi*real(na*nb*nc, kind=dp)*c9/(6._dp*rcut**3*cell%deth**2)
2072  END DO
2073  END IF
2074  END DO
2075  END DO
2076  IF (use_virial) THEN
2077  IF (para_env%is_source()) THEN
2078  DO i = 1, 3
2079  virial%pv_virial(i, i) = virial%pv_virial(i, i) + (elrc6 + elrc8 + 2._dp*elrc9)
2080  END DO
2081  END IF
2082  END IF
2083  DEALLOCATE (cnkind)
2084  END IF
2085  END IF
2086 
2087  IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
2088  dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
2089  DEALLOCATE (cnumbers)
2090  IF (dispersion_env%doabc .AND. dispersion_env%c9cnst) THEN
2091  DEALLOCATE (cnumfix)
2092  END IF
2093  IF (calculate_forces .OR. debugall) THEN
2094  DO iatom = 1, natom
2095  DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
2096  END DO
2097  DEALLOCATE (dcnum)
2098  ELSE
2099  DEALLOCATE (dcnum)
2100  END IF
2101  END IF
2102 
2103  ! set dispersion energy
2104  CALL para_env%sum(evdw)
2105  evdw = evdw + (elrc6 + elrc8 + elrc9)
2106  energy = evdw
2107 
2108  IF ((dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
2109  dispersion_env%pp_type == vdw_pairpot_dftd3bj) .AND. debugall) THEN
2110  CALL para_env%sum(e6tot)
2111  CALL para_env%sum(e8tot)
2112  CALL para_env%sum(e9tot)
2113  CALL para_env%sum(nab)
2114  CALL para_env%sum(nabc)
2115  e6tot = e6tot + elrc6
2116  e8tot = e8tot + elrc8
2117  e9tot = e9tot + elrc9
2118  IF (unit_nr > 0) THEN
2119  WRITE (unit_nr, "(A,F20.0)") " E6 vdW terms :", nab
2120  WRITE (unit_nr, *) " E6 vdW energy [au/kcal] :", e6tot, e6tot*kcalmol
2121  WRITE (unit_nr, *) " E8 vdW energy [au/kcal] :", e8tot, e8tot*kcalmol
2122  WRITE (unit_nr, *) " %E8 on total vdW energy :", e8tot/evdw*100._dp
2123  WRITE (unit_nr, "(A,F20.0)") " E9 vdW terms :", nabc
2124  WRITE (unit_nr, *) " E9 vdW energy [au/kcal] :", e9tot, e9tot*kcalmol
2125  WRITE (unit_nr, *) " %E9 on total vdW energy :", e9tot/evdw*100._dp
2126  IF (dispersion_env%lrc) THEN
2127  WRITE (unit_nr, *) " E LRC C6 [au/kcal] :", elrc6, elrc6*kcalmol
2128  WRITE (unit_nr, *) " E LRC C8 [au/kcal] :", elrc8, elrc8*kcalmol
2129  WRITE (unit_nr, *) " E LRC C9 [au/kcal] :", elrc9, elrc9*kcalmol
2130  END IF
2131  END IF
2132  END IF
2133  IF (unit_nr > 0) THEN
2134  WRITE (unit_nr, *) " Total vdW energy [au] :", evdw
2135  WRITE (unit_nr, *) " Total vdW energy [kcal] :", evdw*kcalmol
2136  WRITE (unit_nr, *)
2137  END IF
2138  IF (calculate_forces .AND. debugall) THEN
2139  IF (unit_nr > 0) THEN
2140  WRITE (unit_nr, *) " Dispersion Forces "
2141  WRITE (unit_nr, *) " Atom Kind Forces "
2142  END IF
2143  gnorm = 0._dp
2144  DO iatom = 1, natom
2145  ikind = kind_of(iatom)
2146  atom_a = atom_of_kind(iatom)
2147  fdij(1:3) = force(ikind)%dispersion(:, atom_a)
2148  CALL para_env%sum(fdij)
2149  gnorm = gnorm + sum(abs(fdij))
2150  IF (unit_nr > 0) WRITE (unit_nr, "(i5,i7,3F20.14)") iatom, ikind, fdij
2151  END DO
2152  IF (unit_nr > 0) THEN
2153  WRITE (unit_nr, *)
2154  WRITE (unit_nr, *) "|G| = ", gnorm
2155  WRITE (unit_nr, *)
2156  END IF
2157  IF (use_virial) THEN
2158  dvirial = virial%pv_virial - dvirial
2159  CALL para_env%sum(dvirial)
2160  IF (unit_nr > 0) THEN
2161  WRITE (unit_nr, *) "Stress Tensor (dispersion)"
2162  WRITE (unit_nr, "(3G20.12)") dvirial
2163  WRITE (unit_nr, *) " Tr(P)/3 : ", (dvirial(1, 1) + dvirial(2, 2) + dvirial(3, 3))/3._dp
2164  WRITE (unit_nr, *)
2165  END IF
2166  END IF
2167  END IF
2168 
2169  DEALLOCATE (dodisp, ghost, floating, atomnumber, rcpbc, radd2, c6d2)
2170 
2171  IF (domol) THEN
2172  DEALLOCATE (atom2mol)
2173  END IF
2174 
2175  IF (calculate_forces .AND. use_virial) THEN
2176  virial%pv_vdw = virial%pv_vdw + (virial%pv_virial - pv_loc)
2177  END IF
2178 
2179  IF (ASSOCIATED(dispersion_env%dftd_section)) THEN
2180  CALL cp_print_key_finished_output(unit_nr, logger, dispersion_env%dftd_section, "PRINT_DFTD")
2181  END IF
2182 
2183  CALL timestop(handle)
2184 
2185  END SUBROUTINE calculate_dispersion_pairpot
2186 
2187 ! **************************************************************************************************
2188 !> \brief ...
2189 !> \param cell ...
2190 !> \param ncell ...
2191 !> \return ...
2192 ! **************************************************************************************************
2193  FUNCTION cellhash(cell, ncell) RESULT(hash)
2194  INTEGER, DIMENSION(3), INTENT(IN) :: cell, ncell
2195  INTEGER :: hash
2196 
2197  INTEGER :: ix, iy, iz, nx, ny, nz
2198 
2199  cpassert(all(abs(cell) <= ncell))
2200 
2201  ix = cell(1)
2202  IF (ix /= 0) THEN
2203  ix = 2*abs(ix) - (1 + sign(1, ix))/2
2204  END IF
2205  iy = cell(2)
2206  IF (iy /= 0) THEN
2207  iy = 2*abs(iy) - (1 + sign(1, iy))/2
2208  END IF
2209  iz = cell(3)
2210  IF (iz /= 0) THEN
2211  iz = 2*abs(iz) - (1 + sign(1, iz))/2
2212  END IF
2213 
2214  nx = 2*ncell(1) + 1
2215  ny = 2*ncell(2) + 1
2216  nz = 2*ncell(3) + 1
2217 
2218  hash = ix*ny*nz + iy*nz + iz + 1
2219 
2220  END FUNCTION cellhash
2221 ! **************************************************************************************************
2222 !> \brief ...
2223 !> \param z ...
2224 !> \param c6 ...
2225 !> \param r ...
2226 !> \param found ...
2227 ! **************************************************************************************************
2228  SUBROUTINE dftd2_param(z, c6, r, found)
2229 
2230  INTEGER, INTENT(in) :: z
2231  REAL(kind=dp), INTENT(inout) :: c6, r
2232  LOGICAL, INTENT(inout) :: found
2233 
2234  REAL(kind=dp), DIMENSION(54), PARAMETER :: c6val = (/0.14_dp, 0.08_dp, 1.61_dp, 1.61_dp, &
2235  3.13_dp, 1.75_dp, 1.23_dp, 0.70_dp, 0.75_dp, 0.63_dp, 5.71_dp, 5.71_dp, 10.79_dp, 9.23_dp,&
2236  7.84_dp, 5.57_dp, 5.07_dp, 4.61_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, &
2237  10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 16.99_dp, 17.10_dp, &
2238  16.37_dp, 12.64_dp, 12.47_dp, 12.01_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, &
2239  24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 37.32_dp, 38.71_dp, &
2240  38.44_dp, 31.74_dp, 31.50_dp, 29.99_dp/)
2241  REAL(kind=dp), DIMENSION(54), PARAMETER :: rval = (/1.001_dp, 1.012_dp, 0.825_dp, 1.408_dp, &
2242  1.485_dp, 1.452_dp, 1.397_dp, 1.342_dp, 1.287_dp, 1.243_dp, 1.144_dp, 1.364_dp, 1.639_dp, &
2243  1.716_dp, 1.705_dp, 1.683_dp, 1.639_dp, 1.595_dp, 1.485_dp, 1.474_dp, 1.562_dp, 1.562_dp, &
2244  1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.650_dp, &
2245  1.727_dp, 1.760_dp, 1.771_dp, 1.749_dp, 1.727_dp, 1.628_dp, 1.606_dp, 1.639_dp, 1.639_dp, &
2246  1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.672_dp, &
2247  1.804_dp, 1.881_dp, 1.892_dp, 1.892_dp, 1.881_dp/)
2248 
2249 !
2250 ! GRIMME DISPERSION PARAMETERS
2251 ! Stefan Grimme, Semiempirical GGA-Type Density Functional Constructed
2252 ! with a Long-Range Dispersion Correction, J. Comp. Chem. 27: 1787-1799 (2006)
2253 ! doi:10.1002/jcc.20495
2254 !
2255 ! Conversion factor [Jnm^6mol^-1] -> [a.u.] : 17.34527758021901
2256 ! Conversion factor [A] -> [a.u.] : 1.889726132885643
2257 !
2258 ! C6 values in [Jnm^6/mol]
2259 ! vdW radii [A]
2260 
2261  IF (z > 0 .AND. z <= 54) THEN
2262  found = .true.
2263  c6 = c6val(z)*1000._dp*bohr**6/kjmol
2264  r = rval(z)*bohr
2265  ELSE
2266  found = .false.
2267  END IF
2268 
2269  END SUBROUTINE dftd2_param
2270 ! **************************************************************************************************
2271 !> \brief ...
2272 !> \param c6ab ...
2273 !> \param maxci ...
2274 !> \param filename ...
2275 !> \param para_env ...
2276 ! **************************************************************************************************
2277  SUBROUTINE dftd3_c6_param(c6ab, maxci, filename, para_env)
2278 
2279  REAL(kind=dp), DIMENSION(:, :, :, :, :) :: c6ab
2280  INTEGER, DIMENSION(:) :: maxci
2281  CHARACTER(LEN=*) :: filename
2282  TYPE(mp_para_env_type), POINTER :: para_env
2283 
2284  INTEGER :: funit, iadr, iat, jadr, jat, kk, nl, &
2285  nlines, nn
2286  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: pars
2287 
2288  IF (para_env%is_source()) THEN
2289  ! Read the DFT-D3 C6AB parameters from file "filename"
2290  CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
2291  READ (funit, *) nl, nlines
2292  END IF
2293  CALL para_env%bcast(nl)
2294  CALL para_env%bcast(nlines)
2295  ALLOCATE (pars(nl))
2296  IF (para_env%is_source()) THEN
2297  READ (funit, *) pars(1:nl)
2298  CALL close_file(unit_number=funit)
2299  END IF
2300  CALL para_env%bcast(pars)
2301 
2302  ! Store C6AB coefficients in an array
2303  c6ab = -1._dp
2304  maxci = 0
2305  kk = 1
2306  DO nn = 1, nlines
2307  iat = nint(pars(kk + 1))
2308  jat = nint(pars(kk + 2))
2309  CALL limit(iat, jat, iadr, jadr)
2310  maxci(iat) = max(maxci(iat), iadr)
2311  maxci(jat) = max(maxci(jat), jadr)
2312  c6ab(iat, jat, iadr, jadr, 1) = pars(kk)
2313  c6ab(iat, jat, iadr, jadr, 2) = pars(kk + 3)
2314  c6ab(iat, jat, iadr, jadr, 3) = pars(kk + 4)
2315 
2316  c6ab(jat, iat, jadr, iadr, 1) = pars(kk)
2317  c6ab(jat, iat, jadr, iadr, 2) = pars(kk + 4)
2318  c6ab(jat, iat, jadr, iadr, 3) = pars(kk + 3)
2319  kk = (nn*5) + 1
2320  END DO
2321 
2322  DEALLOCATE (pars)
2323 
2324  END SUBROUTINE dftd3_c6_param
2325 
2326 ! **************************************************************************************************
2327 !> \brief ...
2328 !> \param iat ...
2329 !> \param jat ...
2330 !> \param iadr ...
2331 !> \param jadr ...
2332 ! **************************************************************************************************
2333  SUBROUTINE limit(iat, jat, iadr, jadr)
2334  INTEGER :: iat, jat, iadr, jadr
2335 
2336  INTEGER :: i
2337 
2338  iadr = 1
2339  jadr = 1
2340  i = 100
2341  DO WHILE (iat .GT. 100)
2342  iat = iat - 100
2343  iadr = iadr + 1
2344  END DO
2345 
2346  i = 100
2347  DO WHILE (jat .GT. 100)
2348  jat = jat - 100
2349  jadr = jadr + 1
2350  END DO
2351  END SUBROUTINE limit
2352 
2353 ! **************************************************************************************************
2354 !> \brief ...
2355 !> \param rout ...
2356 !> \param rcov ...
2357 !> \param r2r4 ...
2358 ! **************************************************************************************************
2359  SUBROUTINE setr0ab(rout, rcov, r2r4)
2360  ! set cut-off radii
2361  REAL(kind=dp), DIMENSION(:, :) :: rout
2362  REAL(kind=dp), DIMENSION(:) :: rcov, r2r4
2363 
2364  INTEGER :: i, j, k
2365  REAL(kind=dp), DIMENSION(4465) :: r0ab(4465)
2366 
2367  r0ab(1:70) = (/ &
2368  2.1823, 1.8547, 1.7347, 2.9086, 2.5732, 3.4956, 2.3550 &
2369  , 2.5095, 2.9802, 3.0982, 2.5141, 2.3917, 2.9977, 2.9484 &
2370  , 3.2160, 2.4492, 2.2527, 3.1933, 3.0214, 2.9531, 2.9103 &
2371  , 2.3667, 2.1328, 2.8784, 2.7660, 2.7776, 2.7063, 2.6225 &
2372  , 2.1768, 2.0625, 2.6395, 2.6648, 2.6482, 2.5697, 2.4846 &
2373  , 2.4817, 2.0646, 1.9891, 2.5086, 2.6908, 2.6233, 2.4770 &
2374  , 2.3885, 2.3511, 2.2996, 1.9892, 1.9251, 2.4190, 2.5473 &
2375  , 2.4994, 2.4091, 2.3176, 2.2571, 2.1946, 2.1374, 2.9898 &
2376  , 2.6397, 3.6031, 3.1219, 3.7620, 3.2485, 2.9357, 2.7093 &
2377  , 2.5781, 2.4839, 3.7082, 2.5129, 2.7321, 3.1052, 3.2962 &
2378  /)
2379  r0ab(71:140) = (/ &
2380  3.1331, 3.2000, 2.9586, 3.0822, 2.8582, 2.7120, 3.2570 &
2381  , 3.4839, 2.8766, 2.7427, 3.2776, 3.2363, 3.5929, 3.2826 &
2382  , 3.0911, 2.9369, 2.9030, 2.7789, 3.3921, 3.3970, 4.0106 &
2383  , 2.8884, 2.6605, 3.7513, 3.1613, 3.3605, 3.3325, 3.0991 &
2384  , 2.9297, 2.8674, 2.7571, 3.8129, 3.3266, 3.7105, 3.7917 &
2385  , 2.8304, 2.5538, 3.3932, 3.1193, 3.1866, 3.1245, 3.0465 &
2386  , 2.8727, 2.7664, 2.6926, 3.4608, 3.2984, 3.5142, 3.5418 &
2387  , 3.5017, 2.6190, 2.4797, 3.1331, 3.0540, 3.0651, 2.9879 &
2388  , 2.9054, 2.8805, 2.7330, 2.6331, 3.2096, 3.5668, 3.3684 &
2389  , 3.3686, 3.3180, 3.3107, 2.4757, 2.4019, 2.9789, 3.1468 &
2390  /)
2391  r0ab(141:210) = (/ &
2392  2.9768, 2.8848, 2.7952, 2.7457, 2.6881, 2.5728, 3.0574 &
2393  , 3.3264, 3.3562, 3.2529, 3.1916, 3.1523, 3.1046, 2.3725 &
2394  , 2.3289, 2.8760, 2.9804, 2.9093, 2.8040, 2.7071, 2.6386 &
2395  , 2.5720, 2.5139, 2.9517, 3.1606, 3.2085, 3.1692, 3.0982 &
2396  , 3.0352, 2.9730, 2.9148, 3.2147, 2.8315, 3.8724, 3.4621 &
2397  , 3.8823, 3.3760, 3.0746, 2.8817, 2.7552, 2.6605, 3.9740 &
2398  , 3.6192, 3.6569, 3.9586, 3.6188, 3.3917, 3.2479, 3.1434 &
2399  , 4.2411, 2.7597, 3.0588, 3.3474, 3.6214, 3.4353, 3.4729 &
2400  , 3.2487, 3.3200, 3.0914, 2.9403, 3.4972, 3.7993, 3.6773 &
2401  , 3.8678, 3.5808, 3.8243, 3.5826, 3.4156, 3.8765, 4.1035 &
2402  /)
2403  r0ab(211:280) = (/ &
2404  2.7361, 2.9765, 3.2475, 3.5004, 3.4185, 3.4378, 3.2084 &
2405  , 3.2787, 3.0604, 2.9187, 3.4037, 3.6759, 3.6586, 3.8327 &
2406  , 3.5372, 3.7665, 3.5310, 3.3700, 3.7788, 3.9804, 3.8903 &
2407  , 2.6832, 2.9060, 3.2613, 3.4359, 3.3538, 3.3860, 3.1550 &
2408  , 3.2300, 3.0133, 2.8736, 3.4024, 3.6142, 3.5979, 3.5295 &
2409  , 3.4834, 3.7140, 3.4782, 3.3170, 3.7434, 3.9623, 3.8181 &
2410  , 3.7642, 2.6379, 2.8494, 3.1840, 3.4225, 3.2771, 3.3401 &
2411  , 3.1072, 3.1885, 2.9714, 2.8319, 3.3315, 3.5979, 3.5256 &
2412  , 3.4980, 3.4376, 3.6714, 3.4346, 3.2723, 3.6859, 3.8985 &
2413  , 3.7918, 3.7372, 3.7211, 2.9230, 2.6223, 3.4161, 2.8999 &
2414  /)
2415  r0ab(281:350) = (/ &
2416  3.0557, 3.3308, 3.0555, 2.8508, 2.7385, 2.6640, 3.5263 &
2417  , 3.0277, 3.2990, 3.7721, 3.5017, 3.2751, 3.1368, 3.0435 &
2418  , 3.7873, 3.2858, 3.2140, 3.1727, 3.2178, 3.4414, 2.5490 &
2419  , 2.7623, 3.0991, 3.3252, 3.1836, 3.2428, 3.0259, 3.1225 &
2420  , 2.9032, 2.7621, 3.2490, 3.5110, 3.4429, 3.3845, 3.3574 &
2421  , 3.6045, 3.3658, 3.2013, 3.6110, 3.8241, 3.7090, 3.6496 &
2422  , 3.6333, 3.0896, 3.5462, 2.4926, 2.7136, 3.0693, 3.2699 &
2423  , 3.1272, 3.1893, 2.9658, 3.0972, 2.8778, 2.7358, 3.2206 &
2424  , 3.4566, 3.3896, 3.3257, 3.2946, 3.5693, 3.3312, 3.1670 &
2425  , 3.5805, 3.7711, 3.6536, 3.5927, 3.5775, 3.0411, 3.4885 &
2426  /)
2427  r0ab(351:420) = (/ &
2428  3.4421, 2.4667, 2.6709, 3.0575, 3.2357, 3.0908, 3.1537 &
2429  , 2.9235, 3.0669, 2.8476, 2.7054, 3.2064, 3.4519, 3.3593 &
2430  , 3.2921, 3.2577, 3.2161, 3.2982, 3.1339, 3.5606, 3.7582 &
2431  , 3.6432, 3.5833, 3.5691, 3.0161, 3.4812, 3.4339, 3.4327 &
2432  , 2.4515, 2.6338, 3.0511, 3.2229, 3.0630, 3.1265, 2.8909 &
2433  , 3.0253, 2.8184, 2.6764, 3.1968, 3.4114, 3.3492, 3.2691 &
2434  , 3.2320, 3.1786, 3.2680, 3.1036, 3.5453, 3.7259, 3.6090 &
2435  , 3.5473, 3.5327, 3.0018, 3.4413, 3.3907, 3.3593, 3.3462 &
2436  , 2.4413, 2.6006, 3.0540, 3.1987, 3.0490, 3.1058, 2.8643 &
2437  , 2.9948, 2.7908, 2.6491, 3.1950, 3.3922, 3.3316, 3.2585 &
2438  /)
2439  r0ab(421:490) = (/ &
2440  3.2136, 3.1516, 3.2364, 3.0752, 3.5368, 3.7117, 3.5941 &
2441  , 3.5313, 3.5164, 2.9962, 3.4225, 3.3699, 3.3370, 3.3234 &
2442  , 3.3008, 2.4318, 2.5729, 3.0416, 3.1639, 3.0196, 3.0843 &
2443  , 2.8413, 2.7436, 2.7608, 2.6271, 3.1811, 3.3591, 3.3045 &
2444  , 3.2349, 3.1942, 3.1291, 3.2111, 3.0534, 3.5189, 3.6809 &
2445  , 3.5635, 3.5001, 3.4854, 2.9857, 3.3897, 3.3363, 3.3027 &
2446  , 3.2890, 3.2655, 3.2309, 2.8502, 2.6934, 3.2467, 3.1921 &
2447  , 3.5663, 3.2541, 3.0571, 2.9048, 2.8657, 2.7438, 3.3547 &
2448  , 3.3510, 3.9837, 3.6871, 3.4862, 3.3389, 3.2413, 3.1708 &
2449  , 3.6096, 3.6280, 3.6860, 3.5568, 3.4836, 3.2868, 3.3994 &
2450  /)
2451  r0ab(491:560) = (/ &
2452  3.3476, 3.3170, 3.2950, 3.2874, 3.2606, 3.9579, 2.9226 &
2453  , 2.6838, 3.7867, 3.1732, 3.3872, 3.3643, 3.1267, 2.9541 &
2454  , 2.8505, 2.7781, 3.8475, 3.3336, 3.7359, 3.8266, 3.5733 &
2455  , 3.3959, 3.2775, 3.1915, 3.9878, 3.8816, 3.5810, 3.5364 &
2456  , 3.5060, 3.8097, 3.3925, 3.3348, 3.3019, 3.2796, 3.2662 &
2457  , 3.2464, 3.7136, 3.8619, 2.9140, 2.6271, 3.4771, 3.1774 &
2458  , 3.2560, 3.1970, 3.1207, 2.9406, 2.8322, 2.7571, 3.5455 &
2459  , 3.3514, 3.5837, 3.6177, 3.5816, 3.3902, 3.2604, 3.1652 &
2460  , 3.7037, 3.6283, 3.5858, 3.5330, 3.4884, 3.5789, 3.4094 &
2461  , 3.3473, 3.3118, 3.2876, 3.2707, 3.2521, 3.5570, 3.6496 &
2462  /)
2463  r0ab(561:630) = (/ &
2464  3.6625, 2.7300, 2.5870, 3.2471, 3.1487, 3.1667, 3.0914 &
2465  , 3.0107, 2.9812, 2.8300, 2.7284, 3.3259, 3.3182, 3.4707 &
2466  , 3.4748, 3.4279, 3.4182, 3.2547, 3.1353, 3.5116, 3.9432 &
2467  , 3.8828, 3.8303, 3.7880, 3.3760, 3.7218, 3.3408, 3.3059 &
2468  , 3.2698, 3.2446, 3.2229, 3.4422, 3.5023, 3.5009, 3.5268 &
2469  , 2.6026, 2.5355, 3.1129, 3.2863, 3.1029, 3.0108, 2.9227 &
2470  , 2.8694, 2.8109, 2.6929, 3.1958, 3.4670, 3.4018, 3.3805 &
2471  , 3.3218, 3.2815, 3.2346, 3.0994, 3.3937, 3.7266, 3.6697 &
2472  , 3.6164, 3.5730, 3.2522, 3.5051, 3.4686, 3.4355, 3.4084 &
2473  , 3.3748, 3.3496, 3.3692, 3.4052, 3.3910, 3.3849, 3.3662 &
2474  /)
2475  r0ab(631:700) = (/ &
2476  2.5087, 2.4814, 3.0239, 3.1312, 3.0535, 2.9457, 2.8496 &
2477  , 2.7780, 2.7828, 2.6532, 3.1063, 3.3143, 3.3549, 3.3120 &
2478  , 3.2421, 3.1787, 3.1176, 3.0613, 3.3082, 3.5755, 3.5222 &
2479  , 3.4678, 3.4231, 3.1684, 3.3528, 3.3162, 3.2827, 3.2527 &
2480  , 3.2308, 3.2029, 3.3173, 3.3343, 3.3092, 3.2795, 3.2452 &
2481  , 3.2096, 3.2893, 2.8991, 4.0388, 3.6100, 3.9388, 3.4475 &
2482  , 3.1590, 2.9812, 2.8586, 2.7683, 4.1428, 3.7911, 3.8225 &
2483  , 4.0372, 3.7059, 3.4935, 3.3529, 3.2492, 4.4352, 4.0826 &
2484  , 3.9733, 3.9254, 3.8646, 3.9315, 3.7837, 3.7465, 3.7211 &
2485  , 3.7012, 3.6893, 3.6676, 3.7736, 4.0660, 3.7926, 3.6158 &
2486  /)
2487  r0ab(701:770) = (/ &
2488  3.5017, 3.4166, 4.6176, 2.8786, 3.1658, 3.5823, 3.7689 &
2489  , 3.5762, 3.5789, 3.3552, 3.4004, 3.1722, 3.0212, 3.7241 &
2490  , 3.9604, 3.8500, 3.9844, 3.7035, 3.9161, 3.6751, 3.5075 &
2491  , 4.1151, 4.2877, 4.1579, 4.1247, 4.0617, 3.4874, 3.9848 &
2492  , 3.9280, 3.9079, 3.8751, 3.8604, 3.8277, 3.8002, 3.9981 &
2493  , 3.7544, 4.0371, 3.8225, 3.6718, 4.3092, 4.4764, 2.8997 &
2494  , 3.0953, 3.4524, 3.6107, 3.6062, 3.5783, 3.3463, 3.3855 &
2495  , 3.1746, 3.0381, 3.6019, 3.7938, 3.8697, 3.9781, 3.6877 &
2496  , 3.8736, 3.6451, 3.4890, 3.9858, 4.1179, 4.0430, 3.9563 &
2497  , 3.9182, 3.4002, 3.8310, 3.7716, 3.7543, 3.7203, 3.7053 &
2498  /)
2499  r0ab(771:840) = (/ &
2500  3.6742, 3.8318, 3.7631, 3.7392, 3.9892, 3.7832, 3.6406 &
2501  , 4.1701, 4.3016, 4.2196, 2.8535, 3.0167, 3.3978, 3.5363 &
2502  , 3.5393, 3.5301, 3.2960, 3.3352, 3.1287, 2.9967, 3.6659 &
2503  , 3.7239, 3.8070, 3.7165, 3.6368, 3.8162, 3.5885, 3.4336 &
2504  , 3.9829, 4.0529, 3.9584, 3.9025, 3.8607, 3.3673, 3.7658 &
2505  , 3.7035, 3.6866, 3.6504, 3.6339, 3.6024, 3.7708, 3.7283 &
2506  , 3.6896, 3.9315, 3.7250, 3.5819, 4.1457, 4.2280, 4.1130 &
2507  , 4.0597, 3.0905, 2.7998, 3.6448, 3.0739, 3.2996, 3.5262 &
2508  , 3.2559, 3.0518, 2.9394, 2.8658, 3.7514, 3.2295, 3.5643 &
2509  , 3.7808, 3.6931, 3.4723, 3.3357, 3.2429, 4.0280, 3.5589 &
2510  /)
2511  r0ab(841:910) = (/ &
2512  3.4636, 3.4994, 3.4309, 3.6177, 3.2946, 3.2376, 3.2050 &
2513  , 3.1847, 3.1715, 3.1599, 3.5555, 3.8111, 3.7693, 3.5718 &
2514  , 3.4498, 3.3662, 4.1608, 3.7417, 3.6536, 3.6154, 3.8596 &
2515  , 3.0301, 2.7312, 3.5821, 3.0473, 3.2137, 3.4679, 3.1975 &
2516  , 2.9969, 2.8847, 2.8110, 3.6931, 3.2076, 3.4943, 3.5956 &
2517  , 3.6379, 3.4190, 3.2808, 3.1860, 3.9850, 3.5105, 3.4330 &
2518  , 3.3797, 3.4155, 3.6033, 3.2737, 3.2145, 3.1807, 3.1596 &
2519  , 3.1461, 3.1337, 3.4812, 3.6251, 3.7152, 3.5201, 3.3966 &
2520  , 3.3107, 4.1128, 3.6899, 3.6082, 3.5604, 3.7834, 3.7543 &
2521  , 2.9189, 2.6777, 3.4925, 2.9648, 3.1216, 3.2940, 3.0975 &
2522  /)
2523  r0ab(911:980) = (/ &
2524  2.9757, 2.8493, 2.7638, 3.6085, 3.1214, 3.4006, 3.4793 &
2525  , 3.5147, 3.3806, 3.2356, 3.1335, 3.9144, 3.4183, 3.3369 &
2526  , 3.2803, 3.2679, 3.4871, 3.1714, 3.1521, 3.1101, 3.0843 &
2527  , 3.0670, 3.0539, 3.3890, 3.5086, 3.5895, 3.4783, 3.3484 &
2528  , 3.2559, 4.0422, 3.5967, 3.5113, 3.4576, 3.6594, 3.6313 &
2529  , 3.5690, 2.8578, 2.6334, 3.4673, 2.9245, 3.0732, 3.2435 &
2530  , 3.0338, 2.9462, 2.8143, 2.7240, 3.5832, 3.0789, 3.3617 &
2531  , 3.4246, 3.4505, 3.3443, 3.1964, 3.0913, 3.8921, 3.3713 &
2532  , 3.2873, 3.2281, 3.2165, 3.4386, 3.1164, 3.1220, 3.0761 &
2533  , 3.0480, 3.0295, 3.0155, 3.3495, 3.4543, 3.5260, 3.4413 &
2534  /)
2535  r0ab(981:1050) = (/ &
2536  3.3085, 3.2134, 4.0170, 3.5464, 3.4587, 3.4006, 3.6027 &
2537  , 3.5730, 3.4945, 3.4623, 2.8240, 2.5960, 3.4635, 2.9032 &
2538  , 3.0431, 3.2115, 2.9892, 2.9148, 2.7801, 2.6873, 3.5776 &
2539  , 3.0568, 3.3433, 3.3949, 3.4132, 3.3116, 3.1616, 3.0548 &
2540  , 3.8859, 3.3719, 3.2917, 3.2345, 3.2274, 3.4171, 3.1293 &
2541  , 3.0567, 3.0565, 3.0274, 3.0087, 2.9939, 3.3293, 3.4249 &
2542  , 3.4902, 3.4091, 3.2744, 3.1776, 4.0078, 3.5374, 3.4537 &
2543  , 3.3956, 3.5747, 3.5430, 3.4522, 3.4160, 3.3975, 2.8004 &
2544  , 2.5621, 3.4617, 2.9154, 3.0203, 3.1875, 2.9548, 2.8038 &
2545  , 2.7472, 2.6530, 3.5736, 3.0584, 3.3304, 3.3748, 3.3871 &
2546  /)
2547  r0ab(1051:1120) = (/ &
2548  3.2028, 3.1296, 3.0214, 3.8796, 3.3337, 3.2492, 3.1883 &
2549  , 3.1802, 3.4050, 3.0756, 3.0478, 3.0322, 3.0323, 3.0163 &
2550  , 3.0019, 3.3145, 3.4050, 3.4656, 3.3021, 3.2433, 3.1453 &
2551  , 3.9991, 3.5017, 3.4141, 3.3520, 3.5583, 3.5251, 3.4243 &
2552  , 3.3851, 3.3662, 3.3525, 2.7846, 2.5324, 3.4652, 2.8759 &
2553  , 3.0051, 3.1692, 2.9273, 2.7615, 2.7164, 2.6212, 3.5744 &
2554  , 3.0275, 3.3249, 3.3627, 3.3686, 3.1669, 3.0584, 2.9915 &
2555  , 3.8773, 3.3099, 3.2231, 3.1600, 3.1520, 3.4023, 3.0426 &
2556  , 3.0099, 2.9920, 2.9809, 2.9800, 2.9646, 3.3068, 3.3930 &
2557  , 3.4486, 3.2682, 3.1729, 3.1168, 3.9952, 3.4796, 3.3901 &
2558  /)
2559  r0ab(1121:1190) = (/ &
2560  3.3255, 3.5530, 3.5183, 3.4097, 3.3683, 3.3492, 3.3360 &
2561  , 3.3308, 2.5424, 2.6601, 3.2555, 3.2807, 3.1384, 3.1737 &
2562  , 2.9397, 2.8429, 2.8492, 2.7225, 3.3875, 3.4910, 3.4520 &
2563  , 3.3608, 3.3036, 3.2345, 3.2999, 3.1487, 3.7409, 3.8392 &
2564  , 3.7148, 3.6439, 3.6182, 3.1753, 3.5210, 3.4639, 3.4265 &
2565  , 3.4075, 3.3828, 3.3474, 3.4071, 3.3754, 3.3646, 3.3308 &
2566  , 3.4393, 3.2993, 3.8768, 3.9891, 3.8310, 3.7483, 3.3417 &
2567  , 3.3019, 3.2250, 3.1832, 3.1578, 3.1564, 3.1224, 3.4620 &
2568  , 2.9743, 2.8058, 3.4830, 3.3474, 3.6863, 3.3617, 3.1608 &
2569  , 3.0069, 2.9640, 2.8427, 3.5885, 3.5219, 4.1314, 3.8120 &
2570  /)
2571  r0ab(1191:1260) = (/ &
2572  3.6015, 3.4502, 3.3498, 3.2777, 3.8635, 3.8232, 3.8486 &
2573  , 3.7215, 3.6487, 3.4724, 3.5627, 3.5087, 3.4757, 3.4517 &
2574  , 3.4423, 3.4139, 4.1028, 3.8388, 3.6745, 3.5562, 3.4806 &
2575  , 3.4272, 4.0182, 3.9991, 4.0007, 3.9282, 3.7238, 3.6498 &
2576  , 3.5605, 3.5211, 3.5009, 3.4859, 3.4785, 3.5621, 4.2623 &
2577  , 3.0775, 2.8275, 4.0181, 3.3385, 3.5379, 3.5036, 3.2589 &
2578  , 3.0804, 3.0094, 2.9003, 4.0869, 3.5088, 3.9105, 3.9833 &
2579  , 3.7176, 3.5323, 3.4102, 3.3227, 4.2702, 4.0888, 3.7560 &
2580  , 3.7687, 3.6681, 3.6405, 3.5569, 3.4990, 3.4659, 3.4433 &
2581  , 3.4330, 3.4092, 3.8867, 4.0190, 3.7961, 3.6412, 3.5405 &
2582  /)
2583  r0ab(1261:1330) = (/ &
2584  3.4681, 4.3538, 4.2136, 3.9381, 3.8912, 3.9681, 3.7909 &
2585  , 3.6774, 3.6262, 3.5999, 3.5823, 3.5727, 3.5419, 4.0245 &
2586  , 4.1874, 3.0893, 2.7917, 3.7262, 3.3518, 3.4241, 3.5433 &
2587  , 3.2773, 3.0890, 2.9775, 2.9010, 3.8048, 3.5362, 3.7746 &
2588  , 3.7911, 3.7511, 3.5495, 3.4149, 3.3177, 4.0129, 3.8370 &
2589  , 3.7739, 3.7125, 3.7152, 3.7701, 3.5813, 3.5187, 3.4835 &
2590  , 3.4595, 3.4439, 3.4242, 3.7476, 3.8239, 3.8346, 3.6627 &
2591  , 3.5479, 3.4639, 4.1026, 3.9733, 3.9292, 3.8667, 3.9513 &
2592  , 3.8959, 3.7698, 3.7089, 3.6765, 3.6548, 3.6409, 3.5398 &
2593  , 3.8759, 3.9804, 4.0150, 2.9091, 2.7638, 3.5066, 3.3377 &
2594  /)
2595  r0ab(1331:1400) = (/ &
2596  3.3481, 3.2633, 3.1810, 3.1428, 2.9872, 2.8837, 3.5929 &
2597  , 3.5183, 3.6729, 3.6596, 3.6082, 3.5927, 3.4224, 3.2997 &
2598  , 3.8190, 4.1865, 4.1114, 4.0540, 3.6325, 3.5697, 3.5561 &
2599  , 3.5259, 3.4901, 3.4552, 3.4315, 3.4091, 3.6438, 3.6879 &
2600  , 3.6832, 3.7043, 3.5557, 3.4466, 3.9203, 4.2919, 4.2196 &
2601  , 4.1542, 3.7573, 3.7039, 3.6546, 3.6151, 3.5293, 3.4849 &
2602  , 3.4552, 3.5192, 3.7673, 3.8359, 3.8525, 3.8901, 2.7806 &
2603  , 2.7209, 3.3812, 3.4958, 3.2913, 3.1888, 3.0990, 3.0394 &
2604  , 2.9789, 2.8582, 3.4716, 3.6883, 3.6105, 3.5704, 3.5059 &
2605  , 3.4619, 3.4138, 3.2742, 3.7080, 3.9773, 3.9010, 3.8409 &
2606  /)
2607  r0ab(1401:1470) = (/ &
2608  3.7944, 3.4465, 3.7235, 3.6808, 3.6453, 3.6168, 3.5844 &
2609  , 3.5576, 3.5772, 3.5959, 3.5768, 3.5678, 3.5486, 3.4228 &
2610  , 3.8107, 4.0866, 4.0169, 3.9476, 3.6358, 3.5800, 3.5260 &
2611  , 3.4838, 3.4501, 3.4204, 3.3553, 3.6487, 3.6973, 3.7398 &
2612  , 3.7405, 3.7459, 3.7380, 2.6848, 2.6740, 3.2925, 3.3386 &
2613  , 3.2473, 3.1284, 3.0301, 2.9531, 2.9602, 2.8272, 3.3830 &
2614  , 3.5358, 3.5672, 3.5049, 3.4284, 3.3621, 3.3001, 3.2451 &
2615  , 3.6209, 3.8299, 3.7543, 3.6920, 3.6436, 3.3598, 3.5701 &
2616  , 3.5266, 3.4904, 3.4590, 3.4364, 3.4077, 3.5287, 3.5280 &
2617  , 3.4969, 3.4650, 3.4304, 3.3963, 3.7229, 3.9402, 3.8753 &
2618  /)
2619  r0ab(1471:1540) = (/ &
2620  3.8035, 3.5499, 3.4913, 3.4319, 3.3873, 3.3520, 3.3209 &
2621  , 3.2948, 3.5052, 3.6465, 3.6696, 3.6577, 3.6388, 3.6142 &
2622  , 3.5889, 3.3968, 3.0122, 4.2241, 3.7887, 4.0049, 3.5384 &
2623  , 3.2698, 3.1083, 2.9917, 2.9057, 4.3340, 3.9900, 4.6588 &
2624  , 4.1278, 3.8125, 3.6189, 3.4851, 3.3859, 4.6531, 4.3134 &
2625  , 4.2258, 4.1309, 4.0692, 4.0944, 3.9850, 3.9416, 3.9112 &
2626  , 3.8873, 3.8736, 3.8473, 4.6027, 4.1538, 3.8994, 3.7419 &
2627  , 3.6356, 3.5548, 4.8353, 4.5413, 4.3891, 4.3416, 4.3243 &
2628  , 4.2753, 4.2053, 4.1790, 4.1685, 4.1585, 4.1536, 4.0579 &
2629  , 4.1980, 4.4564, 4.2192, 4.0528, 3.9489, 3.8642, 5.0567 &
2630  /)
2631  r0ab(1541:1610) = (/ &
2632  3.0630, 3.3271, 4.0432, 4.0046, 4.1555, 3.7426, 3.5130 &
2633  , 3.5174, 3.2884, 3.1378, 4.1894, 4.2321, 4.1725, 4.1833 &
2634  , 3.8929, 4.0544, 3.8118, 3.6414, 4.6373, 4.6268, 4.4750 &
2635  , 4.4134, 4.3458, 3.8582, 4.2583, 4.1898, 4.1562, 4.1191 &
2636  , 4.1069, 4.0639, 4.1257, 4.1974, 3.9532, 4.1794, 3.9660 &
2637  , 3.8130, 4.8160, 4.8272, 4.6294, 4.5840, 4.0770, 4.0088 &
2638  , 3.9103, 3.8536, 3.8324, 3.7995, 3.7826, 4.2294, 4.3380 &
2639  , 4.4352, 4.1933, 4.4580, 4.2554, 4.1072, 5.0454, 5.1814 &
2640  , 3.0632, 3.2662, 3.6432, 3.8088, 3.7910, 3.7381, 3.5093 &
2641  , 3.5155, 3.3047, 3.1681, 3.7871, 3.9924, 4.0637, 4.1382 &
2642  /)
2643  r0ab(1611:1680) = (/ &
2644  3.8591, 4.0164, 3.7878, 3.6316, 4.1741, 4.3166, 4.2395 &
2645  , 4.1831, 4.1107, 3.5857, 4.0270, 3.9676, 3.9463, 3.9150 &
2646  , 3.9021, 3.8708, 4.0240, 4.1551, 3.9108, 4.1337, 3.9289 &
2647  , 3.7873, 4.3666, 4.5080, 4.4232, 4.3155, 3.8461, 3.8007 &
2648  , 3.6991, 3.6447, 3.6308, 3.5959, 3.5749, 4.0359, 4.3124 &
2649  , 4.3539, 4.1122, 4.3772, 4.1785, 4.0386, 4.7004, 4.8604 &
2650  , 4.6261, 2.9455, 3.2470, 3.6108, 3.8522, 3.6625, 3.6598 &
2651  , 3.4411, 3.4660, 3.2415, 3.0944, 3.7514, 4.0397, 3.9231 &
2652  , 4.0561, 3.7860, 3.9845, 3.7454, 3.5802, 4.1366, 4.3581 &
2653  , 4.2351, 4.2011, 4.1402, 3.5381, 4.0653, 4.0093, 3.9883 &
2654  /)
2655  r0ab(1681:1750) = (/ &
2656  3.9570, 3.9429, 3.9112, 3.8728, 4.0682, 3.8351, 4.1054 &
2657  , 3.8928, 3.7445, 4.3415, 4.5497, 4.3833, 4.3122, 3.8051 &
2658  , 3.7583, 3.6622, 3.6108, 3.5971, 3.5628, 3.5408, 4.0780 &
2659  , 4.0727, 4.2836, 4.0553, 4.3647, 4.1622, 4.0178, 4.5802 &
2660  , 4.9125, 4.5861, 4.6201, 2.9244, 3.2241, 3.5848, 3.8293 &
2661  , 3.6395, 3.6400, 3.4204, 3.4499, 3.2253, 3.0779, 3.7257 &
2662  , 4.0170, 3.9003, 4.0372, 3.7653, 3.9672, 3.7283, 3.5630 &
2663  , 4.1092, 4.3347, 4.2117, 4.1793, 4.1179, 3.5139, 4.0426 &
2664  , 3.9867, 3.9661, 3.9345, 3.9200, 3.8883, 3.8498, 4.0496 &
2665  , 3.8145, 4.0881, 3.8756, 3.7271, 4.3128, 4.5242, 4.3578 &
2666  /)
2667  r0ab(1751:1820) = (/ &
2668  4.2870, 3.7796, 3.7318, 3.6364, 3.5854, 3.5726, 3.5378 &
2669  , 3.5155, 4.0527, 4.0478, 4.2630, 4.0322, 4.3449, 4.1421 &
2670  , 3.9975, 4.5499, 4.8825, 4.5601, 4.5950, 4.5702, 2.9046 &
2671  , 3.2044, 3.5621, 3.8078, 3.6185, 3.6220, 3.4019, 3.4359 &
2672  , 3.2110, 3.0635, 3.7037, 3.9958, 3.8792, 4.0194, 3.7460 &
2673  , 3.9517, 3.7128, 3.5474, 4.0872, 4.3138, 4.1906, 4.1593 &
2674  , 4.0973, 3.4919, 4.0216, 3.9657, 3.9454, 3.9134, 3.8986 &
2675  , 3.8669, 3.8289, 4.0323, 3.7954, 4.0725, 3.8598, 3.7113 &
2676  , 4.2896, 4.5021, 4.3325, 4.2645, 3.7571, 3.7083, 3.6136 &
2677  , 3.5628, 3.5507, 3.5155, 3.4929, 4.0297, 4.0234, 4.2442 &
2678  /)
2679  r0ab(1821:1890) = (/ &
2680  4.0112, 4.3274, 4.1240, 3.9793, 4.5257, 4.8568, 4.5353 &
2681  , 4.5733, 4.5485, 4.5271, 2.8878, 3.1890, 3.5412, 3.7908 &
2682  , 3.5974, 3.6078, 3.3871, 3.4243, 3.1992, 3.0513, 3.6831 &
2683  , 3.9784, 3.8579, 4.0049, 3.7304, 3.9392, 3.7002, 3.5347 &
2684  , 4.0657, 4.2955, 4.1705, 4.1424, 4.0800, 3.4717, 4.0043 &
2685  , 3.9485, 3.9286, 3.8965, 3.8815, 3.8500, 3.8073, 4.0180 &
2686  , 3.7796, 4.0598, 3.8470, 3.6983, 4.2678, 4.4830, 4.3132 &
2687  , 4.2444, 3.7370, 3.6876, 3.5935, 3.5428, 3.5314, 3.4958 &
2688  , 3.4730, 4.0117, 4.0043, 4.2287, 3.9939, 4.3134, 4.1096 &
2689  , 3.9646, 4.5032, 4.8356, 4.5156, 4.5544, 4.5297, 4.5083 &
2690  /)
2691  r0ab(1891:1960) = (/ &
2692  4.4896, 2.8709, 3.1737, 3.5199, 3.7734, 3.5802, 3.5934 &
2693  , 3.3724, 3.4128, 3.1877, 3.0396, 3.6624, 3.9608, 3.8397 &
2694  , 3.9893, 3.7145, 3.9266, 3.6877, 3.5222, 4.0448, 4.2771 &
2695  , 4.1523, 4.1247, 4.0626, 3.4530, 3.9866, 3.9310, 3.9115 &
2696  , 3.8792, 3.8641, 3.8326, 3.7892, 4.0025, 3.7636, 4.0471 &
2697  , 3.8343, 3.6854, 4.2464, 4.4635, 4.2939, 4.2252, 3.7169 &
2698  , 3.6675, 3.5739, 3.5235, 3.5126, 3.4768, 3.4537, 3.9932 &
2699  , 3.9854, 4.2123, 3.9765, 4.2992, 4.0951, 3.9500, 4.4811 &
2700  , 4.8135, 4.4959, 4.5351, 4.5105, 4.4891, 4.4705, 4.4515 &
2701  , 2.8568, 3.1608, 3.5050, 3.7598, 3.5665, 3.5803, 3.3601 &
2702  /)
2703  r0ab(1961:2030) = (/ &
2704  3.4031, 3.1779, 3.0296, 3.6479, 3.9471, 3.8262, 3.9773 &
2705  , 3.7015, 3.9162, 3.6771, 3.5115, 4.0306, 4.2634, 4.1385 &
2706  , 4.1116, 4.0489, 3.4366, 3.9732, 3.9176, 3.8983, 3.8659 &
2707  , 3.8507, 3.8191, 3.7757, 3.9907, 3.7506, 4.0365, 3.8235 &
2708  , 3.6745, 4.2314, 4.4490, 4.2792, 4.2105, 3.7003, 3.6510 &
2709  , 3.5578, 3.5075, 3.4971, 3.4609, 3.4377, 3.9788, 3.9712 &
2710  , 4.1997, 3.9624, 4.2877, 4.0831, 3.9378, 4.4655, 4.7974 &
2711  , 4.4813, 4.5209, 4.4964, 4.4750, 4.4565, 4.4375, 4.4234 &
2712  , 2.6798, 3.0151, 3.2586, 3.5292, 3.5391, 3.4902, 3.2887 &
2713  , 3.3322, 3.1228, 2.9888, 3.4012, 3.7145, 3.7830, 3.6665 &
2714  /)
2715  r0ab(2031:2100) = (/ &
2716  3.5898, 3.8077, 3.5810, 3.4265, 3.7726, 4.0307, 3.9763 &
2717  , 3.8890, 3.8489, 3.2706, 3.7595, 3.6984, 3.6772, 3.6428 &
2718  , 3.6243, 3.5951, 3.7497, 3.6775, 3.6364, 3.9203, 3.7157 &
2719  , 3.5746, 3.9494, 4.2076, 4.1563, 4.0508, 3.5329, 3.4780 &
2720  , 3.3731, 3.3126, 3.2846, 3.2426, 3.2135, 3.7491, 3.9006 &
2721  , 3.8332, 3.8029, 4.1436, 3.9407, 3.7998, 4.1663, 4.5309 &
2722  , 4.3481, 4.2911, 4.2671, 4.2415, 4.2230, 4.2047, 4.1908 &
2723  , 4.1243, 2.5189, 2.9703, 3.3063, 3.6235, 3.4517, 3.3989 &
2724  , 3.2107, 3.2434, 3.0094, 2.8580, 3.4253, 3.8157, 3.7258 &
2725  , 3.6132, 3.5297, 3.7566, 3.5095, 3.3368, 3.7890, 4.1298 &
2726  /)
2727  r0ab(2101:2170) = (/ &
2728  4.0190, 3.9573, 3.9237, 3.2677, 3.8480, 3.8157, 3.7656 &
2729  , 3.7317, 3.7126, 3.6814, 3.6793, 3.6218, 3.5788, 3.8763 &
2730  , 3.6572, 3.5022, 3.9737, 4.3255, 4.1828, 4.1158, 3.5078 &
2731  , 3.4595, 3.3600, 3.3088, 3.2575, 3.2164, 3.1856, 3.8522 &
2732  , 3.8665, 3.8075, 3.7772, 4.1391, 3.9296, 3.7772, 4.2134 &
2733  , 4.7308, 4.3787, 4.3894, 4.3649, 4.3441, 4.3257, 4.3073 &
2734  , 4.2941, 4.1252, 4.2427, 3.0481, 2.9584, 3.6919, 3.5990 &
2735  , 3.8881, 3.4209, 3.1606, 3.1938, 2.9975, 2.8646, 3.8138 &
2736  , 3.7935, 3.7081, 3.9155, 3.5910, 3.4808, 3.4886, 3.3397 &
2737  , 4.1336, 4.1122, 3.9888, 3.9543, 3.8917, 3.5894, 3.8131 &
2738  /)
2739  r0ab(2171:2240) = (/ &
2740  3.7635, 3.7419, 3.7071, 3.6880, 3.6574, 3.6546, 3.9375 &
2741  , 3.6579, 3.5870, 3.6361, 3.5039, 4.3149, 4.2978, 4.1321 &
2742  , 4.1298, 3.8164, 3.7680, 3.7154, 3.6858, 3.6709, 3.6666 &
2743  , 3.6517, 3.8174, 3.8608, 4.1805, 3.9102, 3.8394, 3.8968 &
2744  , 3.7673, 4.5274, 4.6682, 4.3344, 4.3639, 4.3384, 4.3162 &
2745  , 4.2972, 4.2779, 4.2636, 4.0253, 4.1168, 4.1541, 2.8136 &
2746  , 3.0951, 3.4635, 3.6875, 3.4987, 3.5183, 3.2937, 3.3580 &
2747  , 3.1325, 2.9832, 3.6078, 3.8757, 3.7616, 3.9222, 3.6370 &
2748  , 3.8647, 3.6256, 3.4595, 3.9874, 4.1938, 4.0679, 4.0430 &
2749  , 3.9781, 3.3886, 3.9008, 3.8463, 3.8288, 3.7950, 3.7790 &
2750  /)
2751  r0ab(2241:2310) = (/ &
2752  3.7472, 3.7117, 3.9371, 3.6873, 3.9846, 3.7709, 3.6210 &
2753  , 4.1812, 4.3750, 4.2044, 4.1340, 3.6459, 3.5929, 3.5036 &
2754  , 3.4577, 3.4528, 3.4146, 3.3904, 3.9014, 3.9031, 4.1443 &
2755  , 3.8961, 4.2295, 4.0227, 3.8763, 4.4086, 4.7097, 4.4064 &
2756  , 4.4488, 4.4243, 4.4029, 4.3842, 4.3655, 4.3514, 4.1162 &
2757  , 4.2205, 4.1953, 4.2794, 2.8032, 3.0805, 3.4519, 3.6700 &
2758  , 3.4827, 3.5050, 3.2799, 3.3482, 3.1233, 2.9747, 3.5971 &
2759  , 3.8586, 3.7461, 3.9100, 3.6228, 3.8535, 3.6147, 3.4490 &
2760  , 3.9764, 4.1773, 4.0511, 4.0270, 3.9614, 3.3754, 3.8836 &
2761  , 3.8291, 3.8121, 3.7780, 3.7619, 3.7300, 3.6965, 3.9253 &
2762  /)
2763  r0ab(2311:2380) = (/ &
2764  3.6734, 3.9733, 3.7597, 3.6099, 4.1683, 4.3572, 4.1862 &
2765  , 4.1153, 3.6312, 3.5772, 3.4881, 3.4429, 3.4395, 3.4009 &
2766  , 3.3766, 3.8827, 3.8868, 4.1316, 3.8807, 4.2164, 4.0092 &
2767  , 3.8627, 4.3936, 4.6871, 4.3882, 4.4316, 4.4073, 4.3858 &
2768  , 4.3672, 4.3485, 4.3344, 4.0984, 4.2036, 4.1791, 4.2622 &
2769  , 4.2450, 2.7967, 3.0689, 3.4445, 3.6581, 3.4717, 3.4951 &
2770  , 3.2694, 3.3397, 3.1147, 2.9661, 3.5898, 3.8468, 3.7358 &
2771  , 3.9014, 3.6129, 3.8443, 3.6054, 3.4396, 3.9683, 4.1656 &
2772  , 4.0394, 4.0158, 3.9498, 3.3677, 3.8718, 3.8164, 3.8005 &
2773  , 3.7662, 3.7500, 3.7181, 3.6863, 3.9170, 3.6637, 3.9641 &
2774  /)
2775  r0ab(2381:2450) = (/ &
2776  3.7503, 3.6004, 4.1590, 4.3448, 4.1739, 4.1029, 3.6224 &
2777  , 3.5677, 3.4785, 3.4314, 3.4313, 3.3923, 3.3680, 3.8698 &
2778  , 3.8758, 4.1229, 3.8704, 4.2063, 3.9987, 3.8519, 4.3832 &
2779  , 4.6728, 4.3759, 4.4195, 4.3952, 4.3737, 4.3551, 4.3364 &
2780  , 4.3223, 4.0861, 4.1911, 4.1676, 4.2501, 4.2329, 4.2208 &
2781  , 2.7897, 3.0636, 3.4344, 3.6480, 3.4626, 3.4892, 3.2626 &
2782  , 3.3344, 3.1088, 2.9597, 3.5804, 3.8359, 3.7251, 3.8940 &
2783  , 3.6047, 3.8375, 3.5990, 3.4329, 3.9597, 4.1542, 4.0278 &
2784  , 4.0048, 3.9390, 3.3571, 3.8608, 3.8056, 3.7899, 3.7560 &
2785  , 3.7400, 3.7081, 3.6758, 3.9095, 3.6552, 3.9572, 3.7436 &
2786  /)
2787  r0ab(2451:2520) = (/ &
2788  3.5933, 4.1508, 4.3337, 4.1624, 4.0916, 3.6126, 3.5582 &
2789  , 3.4684, 3.4212, 3.4207, 3.3829, 3.3586, 3.8604, 3.8658 &
2790  , 4.1156, 3.8620, 4.1994, 3.9917, 3.8446, 4.3750, 4.6617 &
2791  , 4.3644, 4.4083, 4.3840, 4.3625, 4.3439, 4.3253, 4.3112 &
2792  , 4.0745, 4.1807, 4.1578, 4.2390, 4.2218, 4.2097, 4.1986 &
2793  , 2.8395, 3.0081, 3.3171, 3.4878, 3.5360, 3.5145, 3.2809 &
2794  , 3.3307, 3.1260, 2.9940, 3.4741, 3.6675, 3.7832, 3.6787 &
2795  , 3.6156, 3.8041, 3.5813, 3.4301, 3.8480, 3.9849, 3.9314 &
2796  , 3.8405, 3.8029, 3.2962, 3.7104, 3.6515, 3.6378, 3.6020 &
2797  , 3.5849, 3.5550, 3.7494, 3.6893, 3.6666, 3.9170, 3.7150 &
2798  /)
2799  r0ab(2521:2590) = (/ &
2800  3.5760, 4.0268, 4.1596, 4.1107, 3.9995, 3.5574, 3.5103 &
2801  , 3.4163, 3.3655, 3.3677, 3.3243, 3.2975, 3.7071, 3.9047 &
2802  , 3.8514, 3.8422, 3.8022, 3.9323, 3.7932, 4.2343, 4.4583 &
2803  , 4.3115, 4.2457, 4.2213, 4.1945, 4.1756, 4.1569, 4.1424 &
2804  , 4.0620, 4.0494, 3.9953, 4.0694, 4.0516, 4.0396, 4.0280 &
2805  , 4.0130, 2.9007, 2.9674, 3.8174, 3.5856, 3.6486, 3.5339 &
2806  , 3.2832, 3.3154, 3.1144, 2.9866, 3.9618, 3.8430, 3.9980 &
2807  , 3.8134, 3.6652, 3.7985, 3.5756, 3.4207, 4.4061, 4.2817 &
2808  , 4.1477, 4.0616, 3.9979, 3.6492, 3.8833, 3.8027, 3.7660 &
2809  , 3.7183, 3.6954, 3.6525, 3.9669, 3.8371, 3.7325, 3.9160 &
2810  /)
2811  r0ab(2591:2660) = (/ &
2812  3.7156, 3.5714, 4.6036, 4.4620, 4.3092, 4.2122, 3.8478 &
2813  , 3.7572, 3.6597, 3.5969, 3.5575, 3.5386, 3.5153, 3.7818 &
2814  , 4.1335, 4.0153, 3.9177, 3.8603, 3.9365, 3.7906, 4.7936 &
2815  , 4.7410, 4.5461, 4.5662, 4.5340, 4.5059, 4.4832, 4.4604 &
2816  , 4.4429, 4.2346, 4.4204, 4.3119, 4.3450, 4.3193, 4.3035 &
2817  , 4.2933, 4.1582, 4.2450, 2.8559, 2.9050, 3.8325, 3.5442 &
2818  , 3.5077, 3.4905, 3.2396, 3.2720, 3.0726, 2.9467, 3.9644 &
2819  , 3.8050, 3.8981, 3.7762, 3.6216, 3.7531, 3.5297, 3.3742 &
2820  , 4.3814, 4.2818, 4.1026, 4.0294, 3.9640, 3.6208, 3.8464 &
2821  , 3.7648, 3.7281, 3.6790, 3.6542, 3.6117, 3.8650, 3.8010 &
2822  /)
2823  r0ab(2661:2730) = (/ &
2824  3.6894, 3.8713, 3.6699, 3.5244, 4.5151, 4.4517, 4.2538 &
2825  , 4.1483, 3.8641, 3.7244, 3.6243, 3.5589, 3.5172, 3.4973 &
2826  , 3.4715, 3.7340, 4.0316, 3.9958, 3.8687, 3.8115, 3.8862 &
2827  , 3.7379, 4.7091, 4.7156, 4.5199, 4.5542, 4.5230, 4.4959 &
2828  , 4.4750, 4.4529, 4.4361, 4.1774, 4.3774, 4.2963, 4.3406 &
2829  , 4.3159, 4.3006, 4.2910, 4.1008, 4.1568, 4.0980, 2.8110 &
2830  , 2.8520, 3.7480, 3.5105, 3.4346, 3.3461, 3.1971, 3.2326 &
2831  , 3.0329, 2.9070, 3.8823, 3.7928, 3.8264, 3.7006, 3.5797 &
2832  , 3.7141, 3.4894, 3.3326, 4.3048, 4.2217, 4.0786, 3.9900 &
2833  , 3.9357, 3.6331, 3.8333, 3.7317, 3.6957, 3.6460, 3.6197 &
2834  /)
2835  r0ab(2731:2800) = (/ &
2836  3.5779, 3.7909, 3.7257, 3.6476, 3.5729, 3.6304, 3.4834 &
2837  , 4.4368, 4.3921, 4.2207, 4.1133, 3.8067, 3.7421, 3.6140 &
2838  , 3.5491, 3.5077, 3.4887, 3.4623, 3.6956, 3.9568, 3.8976 &
2839  , 3.8240, 3.7684, 3.8451, 3.6949, 4.6318, 4.6559, 4.4533 &
2840  , 4.4956, 4.4641, 4.4366, 4.4155, 4.3936, 4.3764, 4.1302 &
2841  , 4.3398, 4.2283, 4.2796, 4.2547, 4.2391, 4.2296, 4.0699 &
2842  , 4.1083, 4.0319, 3.9855, 2.7676, 2.8078, 3.6725, 3.4804 &
2843  , 3.3775, 3.2411, 3.1581, 3.1983, 2.9973, 2.8705, 3.8070 &
2844  , 3.7392, 3.7668, 3.6263, 3.5402, 3.6807, 3.4545, 3.2962 &
2845  , 4.2283, 4.1698, 4.0240, 3.9341, 3.8711, 3.5489, 3.7798 &
2846  /)
2847  r0ab(2801:2870) = (/ &
2848  3.7000, 3.6654, 3.6154, 3.5882, 3.5472, 3.7289, 3.6510 &
2849  , 3.6078, 3.5355, 3.5963, 3.4480, 4.3587, 4.3390, 4.1635 &
2850  , 4.0536, 3.7193, 3.6529, 3.5512, 3.4837, 3.4400, 3.4191 &
2851  , 3.3891, 3.6622, 3.8934, 3.8235, 3.7823, 3.7292, 3.8106 &
2852  , 3.6589, 4.5535, 4.6013, 4.3961, 4.4423, 4.4109, 4.3835 &
2853  , 4.3625, 4.3407, 4.3237, 4.0863, 4.2835, 4.1675, 4.2272 &
2854  , 4.2025, 4.1869, 4.1774, 4.0126, 4.0460, 3.9815, 3.9340 &
2855  , 3.8955, 2.6912, 2.7604, 3.6037, 3.4194, 3.3094, 3.1710 &
2856  , 3.0862, 3.1789, 2.9738, 2.8427, 3.7378, 3.6742, 3.6928 &
2857  , 3.5512, 3.4614, 3.4087, 3.4201, 3.2607, 4.1527, 4.0977 &
2858  /)
2859  r0ab(2871:2940) = (/ &
2860  3.9523, 3.8628, 3.8002, 3.4759, 3.7102, 3.6466, 3.6106 &
2861  , 3.5580, 3.5282, 3.4878, 3.6547, 3.5763, 3.5289, 3.5086 &
2862  , 3.5593, 3.4099, 4.2788, 4.2624, 4.0873, 3.9770, 3.6407 &
2863  , 3.5743, 3.5178, 3.4753, 3.3931, 3.3694, 3.3339, 3.6002 &
2864  , 3.8164, 3.7478, 3.7028, 3.6952, 3.7669, 3.6137, 4.4698 &
2865  , 4.5488, 4.3168, 4.3646, 4.3338, 4.3067, 4.2860, 4.2645 &
2866  , 4.2478, 4.0067, 4.2349, 4.0958, 4.1543, 4.1302, 4.1141 &
2867  , 4.1048, 3.9410, 3.9595, 3.8941, 3.8465, 3.8089, 3.7490 &
2868  , 2.7895, 2.5849, 3.6484, 3.0162, 3.1267, 3.2125, 3.0043 &
2869  , 2.9572, 2.8197, 2.7261, 3.7701, 3.2446, 3.5239, 3.4696 &
2870  /)
2871  r0ab(2941:3010) = (/ &
2872  3.4261, 3.3508, 3.1968, 3.0848, 4.1496, 3.6598, 3.5111 &
2873  , 3.4199, 3.3809, 3.5382, 3.2572, 3.2100, 3.1917, 3.1519 &
2874  , 3.1198, 3.1005, 3.5071, 3.5086, 3.5073, 3.4509, 3.3120 &
2875  , 3.2082, 4.2611, 3.8117, 3.6988, 3.5646, 3.6925, 3.6295 &
2876  , 3.5383, 3.4910, 3.4625, 3.4233, 3.4007, 3.2329, 3.6723 &
2877  , 3.6845, 3.6876, 3.6197, 3.4799, 3.3737, 4.4341, 4.0525 &
2878  , 3.9011, 3.8945, 3.8635, 3.8368, 3.8153, 3.7936, 3.7758 &
2879  , 3.4944, 3.4873, 3.9040, 3.7110, 3.6922, 3.6799, 3.6724 &
2880  , 3.5622, 3.6081, 3.5426, 3.4922, 3.4498, 3.3984, 3.4456 &
2881  , 2.7522, 2.5524, 3.5742, 2.9508, 3.0751, 3.0158, 2.9644 &
2882  /)
2883  r0ab(3011:3080) = (/ &
2884  2.8338, 2.7891, 2.6933, 3.6926, 3.1814, 3.4528, 3.4186 &
2885  , 3.3836, 3.2213, 3.1626, 3.0507, 4.0548, 3.5312, 3.4244 &
2886  , 3.3409, 3.2810, 3.4782, 3.1905, 3.1494, 3.1221, 3.1128 &
2887  , 3.0853, 3.0384, 3.4366, 3.4562, 3.4638, 3.3211, 3.2762 &
2888  , 3.1730, 4.1632, 3.6825, 3.5822, 3.4870, 3.6325, 3.5740 &
2889  , 3.4733, 3.4247, 3.3969, 3.3764, 3.3525, 3.1984, 3.5989 &
2890  , 3.6299, 3.6433, 3.4937, 3.4417, 3.3365, 4.3304, 3.9242 &
2891  , 3.7793, 3.7623, 3.7327, 3.7071, 3.6860, 3.6650, 3.6476 &
2892  , 3.3849, 3.3534, 3.8216, 3.5870, 3.5695, 3.5584, 3.5508 &
2893  , 3.4856, 3.5523, 3.4934, 3.4464, 3.4055, 3.3551, 3.3888 &
2894  /)
2895  r0ab(3081:3150) = (/ &
2896  3.3525, 2.7202, 2.5183, 3.4947, 2.8731, 3.0198, 3.1457 &
2897  , 2.9276, 2.7826, 2.7574, 2.6606, 3.6090, 3.0581, 3.3747 &
2898  , 3.3677, 3.3450, 3.1651, 3.1259, 3.0147, 3.9498, 3.3857 &
2899  , 3.2917, 3.2154, 3.1604, 3.4174, 3.0735, 3.0342, 3.0096 &
2900  , 3.0136, 2.9855, 2.9680, 3.3604, 3.4037, 3.4243, 3.2633 &
2901  , 3.1810, 3.1351, 4.0557, 3.5368, 3.4526, 3.3699, 3.5707 &
2902  , 3.5184, 3.4085, 3.3595, 3.3333, 3.3143, 3.3041, 3.1094 &
2903  , 3.5193, 3.5745, 3.6025, 3.4338, 3.3448, 3.2952, 4.2158 &
2904  , 3.7802, 3.6431, 3.6129, 3.5853, 3.5610, 3.5406, 3.5204 &
2905  , 3.5036, 3.2679, 3.2162, 3.7068, 3.4483, 3.4323, 3.4221 &
2906  /)
2907  r0ab(3151:3220) = (/ &
2908  3.4138, 3.3652, 3.4576, 3.4053, 3.3618, 3.3224, 3.2711 &
2909  , 3.3326, 3.2950, 3.2564, 2.5315, 2.6104, 3.2734, 3.2299 &
2910  , 3.1090, 2.9942, 2.9159, 2.8324, 2.8350, 2.7216, 3.3994 &
2911  , 3.4475, 3.4354, 3.3438, 3.2807, 3.2169, 3.2677, 3.1296 &
2912  , 3.7493, 3.8075, 3.6846, 3.6104, 3.5577, 3.2052, 3.4803 &
2913  , 3.4236, 3.3845, 3.3640, 3.3365, 3.3010, 3.3938, 3.3624 &
2914  , 3.3440, 3.3132, 3.4035, 3.2754, 3.8701, 3.9523, 3.8018 &
2915  , 3.7149, 3.3673, 3.3199, 3.2483, 3.2069, 3.1793, 3.1558 &
2916  , 3.1395, 3.4097, 3.5410, 3.5228, 3.5116, 3.4921, 3.4781 &
2917  , 3.4690, 4.0420, 4.1759, 4.0078, 4.0450, 4.0189, 3.9952 &
2918  /)
2919  r0ab(3221:3290) = (/ &
2920  3.9770, 3.9583, 3.9434, 3.7217, 3.8228, 3.7826, 3.8640 &
2921  , 3.8446, 3.8314, 3.8225, 3.6817, 3.7068, 3.6555, 3.6159 &
2922  , 3.5831, 3.5257, 3.2133, 3.1689, 3.1196, 3.3599, 2.9852 &
2923  , 2.7881, 3.5284, 3.3493, 3.6958, 3.3642, 3.1568, 3.0055 &
2924  , 2.9558, 2.8393, 3.6287, 3.5283, 4.1511, 3.8259, 3.6066 &
2925  , 3.4527, 3.3480, 3.2713, 3.9037, 3.8361, 3.8579, 3.7311 &
2926  , 3.6575, 3.5176, 3.5693, 3.5157, 3.4814, 3.4559, 3.4445 &
2927  , 3.4160, 4.1231, 3.8543, 3.6816, 3.5602, 3.4798, 3.4208 &
2928  , 4.0542, 4.0139, 4.0165, 3.9412, 3.7698, 3.6915, 3.6043 &
2929  , 3.5639, 3.5416, 3.5247, 3.5153, 3.5654, 4.2862, 4.0437 &
2930  /)
2931  r0ab(3291:3360) = (/ &
2932  3.8871, 3.7741, 3.6985, 3.6413, 4.2345, 4.3663, 4.3257 &
2933  , 4.0869, 4.0612, 4.0364, 4.0170, 3.9978, 3.9834, 3.9137 &
2934  , 3.8825, 3.8758, 3.9143, 3.8976, 3.8864, 3.8768, 3.9190 &
2935  , 4.1613, 4.0566, 3.9784, 3.9116, 3.8326, 3.7122, 3.6378 &
2936  , 3.5576, 3.5457, 4.3127, 3.1160, 2.8482, 4.0739, 3.3599 &
2937  , 3.5698, 3.5366, 3.2854, 3.1039, 2.9953, 2.9192, 4.1432 &
2938  , 3.5320, 3.9478, 4.0231, 3.7509, 3.5604, 3.4340, 3.3426 &
2939  , 4.3328, 3.8288, 3.7822, 3.7909, 3.6907, 3.6864, 3.5793 &
2940  , 3.5221, 3.4883, 3.4649, 3.4514, 3.4301, 3.9256, 4.0596 &
2941  , 3.8307, 3.6702, 3.5651, 3.4884, 4.4182, 4.2516, 3.9687 &
2942  /)
2943  r0ab(3361:3430) = (/ &
2944  3.9186, 3.9485, 3.8370, 3.7255, 3.6744, 3.6476, 3.6295 &
2945  , 3.6193, 3.5659, 4.0663, 4.2309, 4.0183, 3.8680, 3.7672 &
2946  , 3.6923, 4.5240, 4.4834, 4.1570, 4.3204, 4.2993, 4.2804 &
2947  , 4.2647, 4.2481, 4.2354, 3.8626, 3.8448, 4.2267, 4.1799 &
2948  , 4.1670, 3.8738, 3.8643, 3.8796, 4.0575, 4.0354, 3.9365 &
2949  , 3.8611, 3.7847, 3.7388, 3.6826, 3.6251, 3.5492, 4.0889 &
2950  , 4.2764, 3.1416, 2.8325, 3.7735, 3.3787, 3.4632, 3.5923 &
2951  , 3.3214, 3.1285, 3.0147, 2.9366, 3.8527, 3.5602, 3.8131 &
2952  , 3.8349, 3.7995, 3.5919, 3.4539, 3.3540, 4.0654, 3.8603 &
2953  , 3.7972, 3.7358, 3.7392, 3.8157, 3.6055, 3.5438, 3.5089 &
2954  /)
2955  r0ab(3431:3500) = (/ &
2956  3.4853, 3.4698, 3.4508, 3.7882, 3.8682, 3.8837, 3.7055 &
2957  , 3.5870, 3.5000, 4.1573, 4.0005, 3.9568, 3.8936, 3.9990 &
2958  , 3.9433, 3.8172, 3.7566, 3.7246, 3.7033, 3.6900, 3.5697 &
2959  , 3.9183, 4.0262, 4.0659, 3.8969, 3.7809, 3.6949, 4.2765 &
2960  , 4.2312, 4.1401, 4.0815, 4.0580, 4.0369, 4.0194, 4.0017 &
2961  , 3.9874, 3.8312, 3.8120, 3.9454, 3.9210, 3.9055, 3.8951 &
2962  , 3.8866, 3.8689, 3.9603, 3.9109, 3.9122, 3.8233, 3.7438 &
2963  , 3.7436, 3.6981, 3.6555, 3.5452, 3.9327, 4.0658, 4.1175 &
2964  , 2.9664, 2.8209, 3.5547, 3.3796, 3.3985, 3.3164, 3.2364 &
2965  , 3.1956, 3.0370, 2.9313, 3.6425, 3.5565, 3.7209, 3.7108 &
2966  /)
2967  r0ab(3501:3570) = (/ &
2968  3.6639, 3.6484, 3.4745, 3.3492, 3.8755, 4.2457, 3.7758 &
2969  , 3.7161, 3.6693, 3.6155, 3.5941, 3.5643, 3.5292, 3.4950 &
2970  , 3.4720, 3.4503, 3.6936, 3.7392, 3.7388, 3.7602, 3.6078 &
2971  , 3.4960, 3.9800, 4.3518, 4.2802, 3.8580, 3.8056, 3.7527 &
2972  , 3.7019, 3.6615, 3.5768, 3.5330, 3.5038, 3.5639, 3.8192 &
2973  , 3.8883, 3.9092, 3.9478, 3.7995, 3.6896, 4.1165, 4.5232 &
2974  , 4.4357, 4.4226, 4.4031, 4.3860, 4.3721, 4.3580, 4.3466 &
2975  , 4.2036, 4.2037, 3.8867, 4.2895, 4.2766, 4.2662, 4.2598 &
2976  , 3.8408, 3.9169, 3.8681, 3.8250, 3.7855, 3.7501, 3.6753 &
2977  , 3.5499, 3.4872, 3.5401, 3.8288, 3.9217, 3.9538, 4.0054 &
2978  /)
2979  r0ab(3571:3640) = (/ &
2980  2.8388, 2.7890, 3.4329, 3.5593, 3.3488, 3.2486, 3.1615 &
2981  , 3.1000, 3.0394, 2.9165, 3.5267, 3.7479, 3.6650, 3.6263 &
2982  , 3.5658, 3.5224, 3.4762, 3.3342, 3.7738, 4.0333, 3.9568 &
2983  , 3.8975, 3.8521, 3.4929, 3.7830, 3.7409, 3.7062, 3.6786 &
2984  , 3.6471, 3.6208, 3.6337, 3.6519, 3.6363, 3.6278, 3.6110 &
2985  , 3.4825, 3.8795, 4.1448, 4.0736, 4.0045, 3.6843, 3.6291 &
2986  , 3.5741, 3.5312, 3.4974, 3.4472, 3.4034, 3.7131, 3.7557 &
2987  , 3.7966, 3.8005, 3.8068, 3.8015, 3.6747, 4.0222, 4.3207 &
2988  , 4.2347, 4.2191, 4.1990, 4.1811, 4.1666, 4.1521, 4.1401 &
2989  , 3.9970, 3.9943, 3.9592, 4.0800, 4.0664, 4.0559, 4.0488 &
2990  /)
2991  r0ab(3641:3710) = (/ &
2992  3.9882, 4.0035, 3.9539, 3.9138, 3.8798, 3.8355, 3.5359 &
2993  , 3.4954, 3.3962, 3.5339, 3.7595, 3.8250, 3.8408, 3.8600 &
2994  , 3.8644, 2.7412, 2.7489, 3.3374, 3.3950, 3.3076, 3.1910 &
2995  , 3.0961, 3.0175, 3.0280, 2.8929, 3.4328, 3.5883, 3.6227 &
2996  , 3.5616, 3.4894, 3.4241, 3.3641, 3.3120, 3.6815, 3.8789 &
2997  , 3.8031, 3.7413, 3.6939, 3.4010, 3.6225, 3.5797, 3.5443 &
2998  , 3.5139, 3.4923, 3.4642, 3.5860, 3.5849, 3.5570, 3.5257 &
2999  , 3.4936, 3.4628, 3.7874, 3.9916, 3.9249, 3.8530, 3.5932 &
3000  , 3.5355, 3.4757, 3.4306, 3.3953, 3.3646, 3.3390, 3.5637 &
3001  , 3.7053, 3.7266, 3.7177, 3.6996, 3.6775, 3.6558, 3.9331 &
3002  /)
3003  r0ab(3711:3780) = (/ &
3004  4.1655, 4.0879, 4.0681, 4.0479, 4.0299, 4.0152, 4.0006 &
3005  , 3.9883, 3.8500, 3.8359, 3.8249, 3.9269, 3.9133, 3.9025 &
3006  , 3.8948, 3.8422, 3.8509, 3.7990, 3.7570, 3.7219, 3.6762 &
3007  , 3.4260, 3.3866, 3.3425, 3.5294, 3.7022, 3.7497, 3.7542 &
3008  , 3.7494, 3.7370, 3.7216, 3.4155, 3.0522, 4.2541, 3.8218 &
3009  , 4.0438, 3.5875, 3.3286, 3.1682, 3.0566, 2.9746, 4.3627 &
3010  , 4.0249, 4.6947, 4.1718, 3.8639, 3.6735, 3.5435, 3.4479 &
3011  , 4.6806, 4.3485, 4.2668, 4.1690, 4.1061, 4.1245, 4.0206 &
3012  , 3.9765, 3.9458, 3.9217, 3.9075, 3.8813, 3.9947, 4.1989 &
3013  , 3.9507, 3.7960, 3.6925, 3.6150, 4.8535, 4.5642, 4.4134 &
3014  /)
3015  r0ab(3781:3850) = (/ &
3016  4.3688, 4.3396, 4.2879, 4.2166, 4.1888, 4.1768, 4.1660 &
3017  , 4.1608, 4.0745, 4.2289, 4.4863, 4.2513, 4.0897, 3.9876 &
3018  , 3.9061, 5.0690, 5.0446, 4.6186, 4.6078, 4.5780, 4.5538 &
3019  , 4.5319, 4.5101, 4.4945, 4.1912, 4.2315, 4.5534, 4.4373 &
3020  , 4.4224, 4.4120, 4.4040, 4.2634, 4.7770, 4.6890, 4.6107 &
3021  , 4.5331, 4.4496, 4.4082, 4.3095, 4.2023, 4.0501, 4.2595 &
3022  , 4.5497, 4.3056, 4.1506, 4.0574, 3.9725, 5.0796, 3.0548 &
3023  , 3.3206, 3.8132, 3.9720, 3.7675, 3.7351, 3.5167, 3.5274 &
3024  , 3.3085, 3.1653, 3.9500, 4.1730, 4.0613, 4.1493, 3.8823 &
3025  , 4.0537, 3.8200, 3.6582, 4.3422, 4.5111, 4.3795, 4.3362 &
3026  /)
3027  r0ab(3851:3920) = (/ &
3028  4.2751, 3.7103, 4.1973, 4.1385, 4.1129, 4.0800, 4.0647 &
3029  , 4.0308, 4.0096, 4.1619, 3.9360, 4.1766, 3.9705, 3.8262 &
3030  , 4.5348, 4.7025, 4.5268, 4.5076, 3.9562, 3.9065, 3.8119 &
3031  , 3.7605, 3.7447, 3.7119, 3.6916, 4.1950, 4.2110, 4.3843 &
3032  , 4.1631, 4.4427, 4.2463, 4.1054, 4.7693, 5.0649, 4.7365 &
3033  , 4.7761, 4.7498, 4.7272, 4.7076, 4.6877, 4.6730, 4.4274 &
3034  , 4.5473, 4.5169, 4.5975, 4.5793, 4.5667, 4.5559, 4.3804 &
3035  , 4.6920, 4.6731, 4.6142, 4.5600, 4.4801, 4.0149, 3.8856 &
3036  , 3.7407, 4.1545, 4.2253, 4.4229, 4.1923, 4.5022, 4.3059 &
3037  , 4.1591, 4.7883, 4.9294, 3.3850, 3.4208, 3.7004, 3.8800 &
3038  /)
3039  r0ab(3921:3990) = (/ &
3040  3.9886, 3.9040, 3.6719, 3.6547, 3.4625, 3.3370, 3.8394 &
3041  , 4.0335, 4.2373, 4.3023, 4.0306, 4.1408, 3.9297, 3.7857 &
3042  , 4.1907, 4.3230, 4.2664, 4.2173, 4.1482, 3.6823, 4.0711 &
3043  , 4.0180, 4.0017, 3.9747, 3.9634, 3.9383, 4.1993, 4.3205 &
3044  , 4.0821, 4.2547, 4.0659, 3.9359, 4.3952, 4.5176, 4.3888 &
3045  , 4.3607, 3.9583, 3.9280, 3.8390, 3.7971, 3.7955, 3.7674 &
3046  , 3.7521, 4.1062, 4.3633, 4.2991, 4.2767, 4.4857, 4.3039 &
3047  , 4.1762, 4.6197, 4.8654, 4.6633, 4.5878, 4.5640, 4.5422 &
3048  , 4.5231, 4.5042, 4.4901, 4.3282, 4.3978, 4.3483, 4.4202 &
3049  , 4.4039, 4.3926, 4.3807, 4.2649, 4.6135, 4.5605, 4.5232 &
3050  /)
3051  r0ab(3991:4060) = (/ &
3052  4.4676, 4.3948, 4.0989, 3.9864, 3.8596, 4.0942, 4.2720 &
3053  , 4.3270, 4.3022, 4.5410, 4.3576, 4.2235, 4.6545, 4.7447 &
3054  , 4.7043, 3.0942, 3.2075, 3.5152, 3.6659, 3.8289, 3.7459 &
3055  , 3.5156, 3.5197, 3.3290, 3.2069, 3.6702, 3.8448, 4.0340 &
3056  , 3.9509, 3.8585, 3.9894, 3.7787, 3.6365, 4.1425, 4.1618 &
3057  , 4.0940, 4.0466, 3.9941, 3.5426, 3.8952, 3.8327, 3.8126 &
3058  , 3.7796, 3.7635, 3.7356, 4.0047, 3.9655, 3.9116, 4.1010 &
3059  , 3.9102, 3.7800, 4.2964, 4.3330, 4.2622, 4.2254, 3.8195 &
3060  , 3.7560, 3.6513, 3.5941, 3.5810, 3.5420, 3.5178, 3.8861 &
3061  , 4.1459, 4.1147, 4.0772, 4.3120, 4.1207, 3.9900, 4.4733 &
3062  /)
3063  r0ab(4061:4130) = (/ &
3064  4.6157, 4.4580, 4.4194, 4.3954, 4.3739, 4.3531, 4.3343 &
3065  , 4.3196, 4.2140, 4.2339, 4.1738, 4.2458, 4.2278, 4.2158 &
3066  , 4.2039, 4.1658, 4.3595, 4.2857, 4.2444, 4.1855, 4.1122 &
3067  , 3.7839, 3.6879, 3.5816, 3.8633, 4.1585, 4.1402, 4.1036 &
3068  , 4.3694, 4.1735, 4.0368, 4.5095, 4.5538, 4.5240, 4.4252 &
3069  , 3.0187, 3.1918, 3.5127, 3.6875, 3.7404, 3.6943, 3.4702 &
3070  , 3.4888, 3.2914, 3.1643, 3.6669, 3.8724, 3.9940, 4.0816 &
3071  , 3.8054, 3.9661, 3.7492, 3.6024, 4.0428, 4.1951, 4.1466 &
3072  , 4.0515, 4.0075, 3.5020, 3.9158, 3.8546, 3.8342, 3.8008 &
3073  , 3.7845, 3.7549, 3.9602, 3.8872, 3.8564, 4.0793, 3.8835 &
3074  /)
3075  r0ab(4131:4200) = (/ &
3076  3.7495, 4.2213, 4.3704, 4.3300, 4.2121, 3.7643, 3.7130 &
3077  , 3.6144, 3.5599, 3.5474, 3.5093, 3.4853, 3.9075, 4.1115 &
3078  , 4.0473, 4.0318, 4.2999, 4.1050, 3.9710, 4.4320, 4.6706 &
3079  , 4.5273, 4.4581, 4.4332, 4.4064, 4.3873, 4.3684, 4.3537 &
3080  , 4.2728, 4.2549, 4.2032, 4.2794, 4.2613, 4.2491, 4.2375 &
3081  , 4.2322, 4.3665, 4.3061, 4.2714, 4.2155, 4.1416, 3.7660 &
3082  , 3.6628, 3.5476, 3.8790, 4.1233, 4.0738, 4.0575, 4.3575 &
3083  , 4.1586, 4.0183, 4.4593, 4.5927, 4.4865, 4.3813, 4.4594 &
3084  , 2.9875, 3.1674, 3.4971, 3.6715, 3.7114, 3.6692, 3.4446 &
3085  , 3.4676, 3.2685, 3.1405, 3.6546, 3.8579, 3.9637, 4.0581 &
3086  /)
3087  r0ab(4201:4270) = (/ &
3088  3.7796, 3.9463, 3.7275, 3.5792, 4.0295, 4.1824, 4.1247 &
3089  , 4.0357, 3.9926, 3.4827, 3.9007, 3.8392, 3.8191, 3.7851 &
3090  , 3.7687, 3.7387, 3.9290, 3.8606, 3.8306, 4.0601, 3.8625 &
3091  , 3.7269, 4.2062, 4.3566, 4.3022, 4.1929, 3.7401, 3.6888 &
3092  , 3.5900, 3.5350, 3.5226, 3.4838, 3.4594, 3.8888, 4.0813 &
3093  , 4.0209, 4.0059, 4.2810, 4.0843, 3.9486, 4.4162, 4.6542 &
3094  , 4.5005, 4.4444, 4.4196, 4.3933, 4.3741, 4.3552, 4.3406 &
3095  , 4.2484, 4.2413, 4.1907, 4.2656, 4.2474, 4.2352, 4.2236 &
3096  , 4.2068, 4.3410, 4.2817, 4.2479, 4.1921, 4.1182, 3.7346 &
3097  , 3.6314, 3.5168, 3.8582, 4.0927, 4.0469, 4.0313, 4.3391 &
3098  /)
3099  r0ab(4271:4340) = (/ &
3100  4.1381, 3.9962, 4.4429, 4.5787, 4.4731, 4.3588, 4.4270 &
3101  , 4.3957, 2.9659, 3.1442, 3.4795, 3.6503, 3.6814, 3.6476 &
3102  , 3.4222, 3.4491, 3.2494, 3.1209, 3.6324, 3.8375, 3.9397 &
3103  , 3.8311, 3.7581, 3.9274, 3.7085, 3.5598, 4.0080, 4.1641 &
3104  , 4.1057, 4.0158, 3.9726, 3.4667, 3.8802, 3.8188, 3.7989 &
3105  , 3.7644, 3.7474, 3.7173, 3.9049, 3.8424, 3.8095, 4.0412 &
3106  , 3.8436, 3.7077, 4.1837, 4.3366, 4.2816, 4.1686, 3.7293 &
3107  , 3.6709, 3.5700, 3.5153, 3.5039, 3.4684, 3.4437, 3.8663 &
3108  , 4.0575, 4.0020, 3.9842, 4.2612, 4.0643, 3.9285, 4.3928 &
3109  , 4.6308, 4.4799, 4.4244, 4.3996, 4.3737, 4.3547, 4.3358 &
3110  /)
3111  r0ab(4341:4410) = (/ &
3112  4.3212, 4.2275, 4.2216, 4.1676, 4.2465, 4.2283, 4.2161 &
3113  , 4.2045, 4.1841, 4.3135, 4.2562, 4.2226, 4.1667, 4.0932 &
3114  , 3.7134, 3.6109, 3.4962, 3.8352, 4.0688, 4.0281, 4.0099 &
3115  , 4.3199, 4.1188, 3.9768, 4.4192, 4.5577, 4.4516, 4.3365 &
3116  , 4.4058, 4.3745, 4.3539, 2.8763, 3.1294, 3.5598, 3.7465 &
3117  , 3.5659, 3.5816, 3.3599, 3.4024, 3.1877, 3.0484, 3.7009 &
3118  , 3.9451, 3.8465, 3.9873, 3.7079, 3.9083, 3.6756, 3.5150 &
3119  , 4.0829, 4.2780, 4.1511, 4.1260, 4.0571, 3.4865, 3.9744 &
3120  , 3.9150, 3.8930, 3.8578, 3.8402, 3.8073, 3.7977, 4.0036 &
3121  , 3.7604, 4.0288, 3.8210, 3.6757, 4.2646, 4.4558, 4.2862 &
3122  /)
3123  r0ab(4411:4465) = (/ &
3124  4.2122, 3.7088, 3.6729, 3.5800, 3.5276, 3.5165, 3.4783 &
3125  , 3.4539, 3.9553, 3.9818, 4.2040, 3.9604, 4.2718, 4.0689 &
3126  , 3.9253, 4.4869, 4.7792, 4.4918, 4.5342, 4.5090, 4.4868 &
3127  , 4.4680, 4.4486, 4.4341, 4.2023, 4.3122, 4.2710, 4.3587 &
3128  , 4.3407, 4.3281, 4.3174, 4.1499, 4.3940, 4.3895, 4.3260 &
3129  , 4.2725, 4.1961, 3.7361, 3.6193, 3.4916, 3.9115, 3.9914 &
3130  , 3.9809, 3.9866, 4.3329, 4.1276, 3.9782, 4.5097, 4.6769 &
3131  , 4.5158, 4.3291, 4.3609, 4.3462, 4.3265, 4.4341 &
3132  /)
3133 
3134  k = 0
3135  DO i = 1, SIZE(rout, 1)
3136  DO j = 1, i
3137  k = k + 1
3138  rout(i, j) = r0ab(k)*bohr
3139  rout(j, i) = r0ab(k)*bohr
3140  END DO
3141  END DO
3142 
3143  ! covalent radii (taken from Pyykko and Atsumi, Chem. Eur. J. 15, 2009, 188-197)
3144  ! values for metals decreased by 10 %
3145  rcov(1:94) = (/ &
3146  0.32, 0.46, 1.20, 0.94, 0.77, 0.75, 0.71, 0.63, 0.64, 0.67 &
3147  , 1.40, 1.25, 1.13, 1.04, 1.10, 1.02, 0.99, 0.96, 1.76, 1.54 &
3148  , 1.33, 1.22, 1.21, 1.10, 1.07, 1.04, 1.00, 0.99, 1.01, 1.09 &
3149  , 1.12, 1.09, 1.15, 1.10, 1.14, 1.17, 1.89, 1.67, 1.47, 1.39 &
3150  , 1.32, 1.24, 1.15, 1.13, 1.13, 1.08, 1.15, 1.23, 1.28, 1.26 &
3151  , 1.26, 1.23, 1.32, 1.31, 2.09, 1.76, 1.62, 1.47, 1.58, 1.57 &
3152  , 1.56, 1.55, 1.51, 1.52, 1.51, 1.50, 1.49, 1.49, 1.48, 1.53 &
3153  , 1.46, 1.37, 1.31, 1.23, 1.18, 1.16, 1.11, 1.12, 1.13, 1.32 &
3154  , 1.30, 1.30, 1.36, 1.31, 1.38, 1.42, 2.01, 1.81, 1.67, 1.58 &
3155  , 1.52, 1.53, 1.54, 1.55 &
3156  /)
3157 
3158  ! PBE0/def2-QZVP atomic values
3159  r2r4(1:94) = (/ &
3160  8.0589, 3.4698, 29.0974, 14.8517, 11.8799, 7.8715, 5.5588, &
3161  4.7566, 3.8025, 3.1036, 26.1552, 17.2304, 17.7210, 12.7442, &
3162  9.5361, 8.1652, 6.7463, 5.6004, 29.2012, 22.3934, 19.0598, &
3163  16.8590, 15.4023, 12.5589, 13.4788, 12.2309, 11.2809, 10.5569, &
3164  10.1428, 9.4907, 13.4606, 10.8544, 8.9386, 8.1350, 7.1251, &
3165  6.1971, 30.0162, 24.4103, 20.3537, 17.4780, 13.5528, 11.8451, &
3166  11.0355, 10.1997, 9.5414, 9.0061, 8.6417, 8.9975, 14.0834, &
3167  11.8333, 10.0179, 9.3844, 8.4110, 7.5152, 32.7622, 27.5708, &
3168  23.1671, 21.6003, 20.9615, 20.4562, 20.1010, 19.7475, 19.4828, &
3169  15.6013, 19.2362, 17.4717, 17.8321, 17.4237, 17.1954, 17.1631, &
3170  14.5716, 15.8758, 13.8989, 12.4834, 11.4421, 10.2671, 8.3549, &
3171  7.8496, 7.3278, 7.4820, 13.5124, 11.6554, 10.0959, 9.7340, &
3172  8.8584, 8.0125, 29.8135, 26.3157, 19.1885, 15.8542, 16.1305, &
3173  15.6161, 15.1226, 16.1576 &
3174  /)
3175 
3176  END SUBROUTINE setr0ab
3177 
3178 ! **************************************************************************************************
3179 !> \brief ...
3180 !> \param cnout ...
3181 ! **************************************************************************************************
3182  SUBROUTINE setcn(cnout)
3183  ! default coordination numbers
3184  REAL(kind=dp), DIMENSION(:) :: cnout
3185 
3186  INTEGER :: n
3187  REAL(kind=dp), DIMENSION(104) :: cn
3188 
3189  cn(1:104) = (/ &
3190  1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 &
3191  , 2.0, 1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 3.0, 2.0, 3.0, 2.0, 2.0, 2.0, 2.0 &
3192  , 3.0, 4.0, 3.0, 4.0, 1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 4.0, 3.0 &
3193  , 2.0, 1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 1.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0 &
3194  , 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 5.0, 6.0, 7.0 &
3195  , 4.0, 4.0, 4.0, 3.0, 2.0, 1.0, 2.0, 3.0, 4.0, 1.0, 0.0, 1.0, 2.0, 3.0, 4.0 &
3196  , 5.0, 6.0, 5.0, 4.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 3.0, 3.0 &
3197  /)
3198 
3199  cnout = 0._dp
3200  n = min(SIZE(cn), SIZE(cnout))
3201  cnout(1:n) = cn(1:n)
3202 
3203  END SUBROUTINE setcn
3204 
3205 ! **************************************************************************************************
3206 !> \brief ...
3207 !> \param rab ...
3208 !> \param rcova ...
3209 !> \param rcovb ...
3210 !> \param k1 ...
3211 !> \param cnab ...
3212 !> \param dcnab ...
3213 ! **************************************************************************************************
3214  SUBROUTINE cnparam_d3(rab, rcova, rcovb, k1, cnab, dcnab)
3215 
3216  REAL(kind=dp), INTENT(IN) :: rab, rcova, rcovb, k1
3217  REAL(kind=dp), INTENT(OUT) :: cnab, dcnab
3218 
3219  REAL(kind=dp) :: dfz, ee, fz, rco, rr
3220 
3221 ! covalent distance in Bohr
3222 
3223  rco = rcova + rcovb
3224  rr = rco/rab
3225  ! counting function exponential has a better long-range behavior
3226  ! than MHGs inverse damping
3227  ! factor k2 already included into rcov
3228  ee = exp(-k1*(rr - 1.0_dp))
3229  ! force the function to zero using a second step function
3230  fz = 0.5_dp*(1.0_dp - tanh(rab - 2.0_dp*rco))
3231  dfz = 0.5_dp*((tanh(rab - 2.0_dp*rco))**2 - 1.0_dp)
3232  cnab = 1.0_dp/(1.0_dp + ee)*fz
3233  dcnab = -cnab*cnab*k1*rr/rab*ee + 1.0_dp/(1.0_dp + ee)*dfz
3234 
3235  END SUBROUTINE cnparam_d3
3236 
3237 ! **************************************************************************************************
3238 !> \brief ...
3239 !> \param rab ...
3240 !> \param rcutab ...
3241 !> \param srn ...
3242 !> \param alpn ...
3243 !> \param rcut ...
3244 !> \param fdab ...
3245 !> \param dfdab ...
3246 ! **************************************************************************************************
3247  SUBROUTINE damping_d3(rab, rcutab, srn, alpn, rcut, fdab, dfdab)
3248 
3249  REAL(kind=dp), INTENT(IN) :: rab, rcutab, srn, alpn, rcut
3250  REAL(kind=dp), INTENT(OUT) :: fdab, dfdab
3251 
3252  REAL(kind=dp) :: a, b, c, d, dd, dfab, dfcc, dz, fab, &
3253  fcc, rl, rr, ru, z, zz
3254 
3255  rl = rcut - 1._dp
3256  ru = rcut
3257  IF (rab >= ru) THEN
3258  fcc = 0._dp
3259  dfcc = 0._dp
3260  ELSEIF (rab <= rl) THEN
3261  fcc = 1._dp
3262  dfcc = 0._dp
3263  ELSE
3264  z = rab*rab - rl*rl
3265  dz = 2._dp*rab
3266  zz = z*z*z
3267  d = ru*ru - rl*rl
3268  dd = d*d*d
3269  a = -10._dp/dd
3270  b = 15._dp/(dd*d)
3271  c = -6._dp/(dd*d*d)
3272  fcc = 1._dp + zz*(a + b*z + c*z*z)
3273  dfcc = zz*dz/z*(3._dp*a + 4._dp*b*z + 5._dp*c*z*z)
3274  END IF
3275 
3276  rr = 6._dp*(rab/(srn*rcutab))**(-alpn)
3277  fab = 1._dp/(1._dp + rr)
3278  dfab = fab*fab*rr*alpn/rab
3279  fdab = fab*fcc
3280  dfdab = dfab*fcc + fab*dfcc
3281 
3282  END SUBROUTINE damping_d3
3283 
3284 ! **************************************************************************************************
3285 !> \brief ...
3286 !> \param maxc ...
3287 !> \param max_elem ...
3288 !> \param c6ab ...
3289 !> \param mxc ...
3290 !> \param iat ...
3291 !> \param jat ...
3292 !> \param nci ...
3293 !> \param ncj ...
3294 !> \param k3 ...
3295 !> \param c6 ...
3296 !> \param dc6a ...
3297 !> \param dc6b ...
3298 ! **************************************************************************************************
3299  SUBROUTINE getc6(maxc, max_elem, c6ab, mxc, iat, jat, nci, ncj, k3, c6, dc6a, dc6b)
3300 
3301  INTEGER, INTENT(IN) :: maxc, max_elem
3302  REAL(kind=dp), INTENT(IN) :: c6ab(max_elem, max_elem, maxc, maxc, 3)
3303  INTEGER, INTENT(IN) :: mxc(max_elem), iat, jat
3304  REAL(kind=dp), INTENT(IN) :: nci, ncj, k3
3305  REAL(kind=dp), INTENT(OUT) :: c6, dc6a, dc6b
3306 
3307  INTEGER :: i, j
3308  REAL(kind=dp) :: c6mem, cn1, cn2, csum, dtmpa, dtmpb, &
3309  dwa, dwb, dza, dzb, r, rsave, rsum, &
3310  tmp1
3311 
3312 ! the exponential is sensitive to numerics
3313 ! when nci or ncj is much larger than cn1/cn2
3314 
3315  c6mem = -1.0e+99_dp
3316  rsave = 1.0e+99_dp
3317  rsum = 0.0_dp
3318  csum = 0.0_dp
3319  dza = 0.0_dp
3320  dzb = 0.0_dp
3321  dwa = 0.0_dp
3322  dwb = 0.0_dp
3323  c6 = 0.0_dp
3324  DO i = 1, mxc(iat)
3325  DO j = 1, mxc(jat)
3326  c6 = c6ab(iat, jat, i, j, 1)
3327  IF (c6 > 0.0_dp) THEN
3328  cn1 = c6ab(iat, jat, i, j, 2)
3329  cn2 = c6ab(iat, jat, i, j, 3)
3330  ! distance
3331  r = (cn1 - nci)**2 + (cn2 - ncj)**2
3332  IF (r < rsave) THEN
3333  rsave = r
3334  c6mem = c6
3335  END IF
3336  tmp1 = exp(k3*r)
3337  dtmpa = -2.0_dp*k3*(cn1 - nci)*tmp1
3338  dtmpb = -2.0_dp*k3*(cn2 - ncj)*tmp1
3339  rsum = rsum + tmp1
3340  csum = csum + tmp1*c6
3341  dza = dza + dtmpa*c6
3342  dwa = dwa + dtmpa
3343  dzb = dzb + dtmpb*c6
3344  dwb = dwb + dtmpb
3345  END IF
3346  END DO
3347  END DO
3348 
3349  IF (c6 == 0.0_dp) c6mem = 0.0_dp
3350 
3351  IF (rsum > 1.0e-66_dp) THEN
3352  c6 = csum/rsum
3353  dc6a = (dza - c6*dwa)/rsum
3354  dc6b = (dzb - c6*dwb)/rsum
3355  ELSE
3356  c6 = c6mem
3357  dc6a = 0._dp
3358  dc6b = 0._dp
3359  END IF
3360 
3361  END SUBROUTINE getc6
3362 
3363 ! **************************************************************************************************
3364 !> \brief ...
3365 !> \param dcnum ...
3366 !> \param para_env ...
3367 ! **************************************************************************************************
3368  SUBROUTINE dcnum_distribute(dcnum, para_env)
3369 
3370  TYPE(dcnum_type), DIMENSION(:) :: dcnum
3371  TYPE(mp_para_env_type), POINTER :: para_env
3372 
3373  INTEGER :: i, ia, ipe, jn, n, natom, ntmax, ntot
3374  INTEGER, ALLOCATABLE, DIMENSION(:) :: list, nloc
3375  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dvals
3376  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rik
3377 
3378 ! we only have to do something if this is a parallel run
3379 
3380  IF (para_env%num_pe > 1) THEN
3381  natom = SIZE(dcnum)
3382  !pack my dcnum data
3383  ALLOCATE (nloc(natom))
3384  DO ia = 1, natom
3385  nloc(ia) = dcnum(ia)%neighbors
3386  END DO
3387  ntot = sum(nloc)
3388  ntmax = ntot
3389  CALL para_env%max(ntmax)
3390  ALLOCATE (list(ntmax), dvals(ntmax), rik(3, ntmax))
3391  list = 0
3392  dvals = 0._dp
3393  rik = 0._dp
3394  i = 0
3395  DO ia = 1, natom
3396  DO jn = 1, nloc(ia)
3397  i = i + 1
3398  list(i) = dcnum(ia)%nlist(jn)
3399  dvals(i) = dcnum(ia)%dvals(jn)
3400  rik(1:3, i) = dcnum(ia)%rik(1:3, jn)
3401  END DO
3402  END DO
3403  DO ipe = 1, para_env%num_pe - 1
3404  !send/receive packed data
3405  CALL para_env%shift(nloc)
3406  !unpack received data
3407  CALL para_env%shift(list)
3408  CALL para_env%shift(dvals)
3409  CALL para_env%shift(rik)
3410  !add data to local dcnum
3411  i = 0
3412  DO ia = 1, natom
3413  n = dcnum(ia)%neighbors + nloc(ia)
3414  IF (n > SIZE(dcnum(ia)%nlist)) THEN
3415  CALL reallocate(dcnum(ia)%nlist, 1, 2*n)
3416  CALL reallocate(dcnum(ia)%dvals, 1, 2*n)
3417  CALL reallocate(dcnum(ia)%rik, 1, 3, 1, 2*n)
3418  END IF
3419  DO jn = 1, nloc(ia)
3420  i = i + 1
3421  n = dcnum(ia)%neighbors + jn
3422  dcnum(ia)%nlist(n) = list(i)
3423  dcnum(ia)%dvals(n) = dvals(i)
3424  dcnum(ia)%rik(1:3, n) = rik(1:3, i)
3425  END DO
3426  dcnum(ia)%neighbors = dcnum(ia)%neighbors + nloc(ia)
3427  END DO
3428  END DO
3429  DEALLOCATE (nloc)
3430  DEALLOCATE (list, dvals, rik)
3431  END IF
3432 
3433  END SUBROUTINE dcnum_distribute
3434 
3435 ! **************************************************************************************************
3436 !> \brief ...
3437 !> \param qs_env ...
3438 !> \param dispersion_env ...
3439 !> \param cnumbers ...
3440 !> \param dcnum ...
3441 !> \param ghost ...
3442 !> \param floating ...
3443 !> \param atomnumber ...
3444 !> \param calculate_forces ...
3445 !> \param debugall ...
3446 ! **************************************************************************************************
3447  SUBROUTINE d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, &
3448  calculate_forces, debugall)
3449 
3450  TYPE(qs_environment_type), POINTER :: qs_env
3451  TYPE(qs_dispersion_type), POINTER :: dispersion_env
3452  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cnumbers
3453  TYPE(dcnum_type), DIMENSION(:), INTENT(INOUT) :: dcnum
3454  LOGICAL, DIMENSION(:), INTENT(IN) :: ghost, floating
3455  INTEGER, DIMENSION(:), INTENT(IN) :: atomnumber
3456  LOGICAL, INTENT(IN) :: calculate_forces, debugall
3457 
3458  CHARACTER(LEN=*), PARAMETER :: routinen = 'd3_cnumber'
3459 
3460  INTEGER :: handle, iatom, ikind, jatom, jkind, &
3461  mepos, natom, ni, nj, num_pe, za, zb
3462  LOGICAL :: exclude_pair
3463  REAL(kind=dp) :: cnab, dcnab, eps_cn, rcc, rcova, rcovb
3464  REAL(kind=dp), DIMENSION(3) :: rij
3465  TYPE(neighbor_list_iterator_p_type), &
3466  DIMENSION(:), POINTER :: nl_iterator
3467  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3468  POINTER :: sab_cn
3469  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3470 
3471 !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
3472 !$ INTEGER :: lock, number_of_locks
3473 
3474  CALL timeset(routinen, handle)
3475 
3476  ! Calculate coordination numbers
3477  NULLIFY (particle_set)
3478  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
3479  natom = SIZE(particle_set)
3480 
3481  eps_cn = dispersion_env%eps_cn
3482  sab_cn => dispersion_env%sab_cn
3483  num_pe = 1
3484 !$ num_pe = omp_get_max_threads()
3485  CALL neighbor_list_iterator_create(nl_iterator, sab_cn, nthread=num_pe)
3486 !$ number_of_locks = natom
3487 
3488 !$OMP PARALLEL DEFAULT( NONE ) &
3489 !$OMP SHARED( locks, number_of_locks, nl_iterator &
3490 !$OMP , ghost, floating, atomnumber, eps_cn &
3491 !$OMP , dispersion_env, calculate_forces, debugall &
3492 !$OMP , dcnum, cnumbers &
3493 !$OMP ) &
3494 !$OMP PRIVATE( mepos &
3495 !$OMP , ikind, jkind, iatom, jatom, rij, rcc &
3496 !$OMP , za, zb, rcova, rcovb, cnab, dcnab &
3497 !$OMP , ni, nj, exclude_pair &
3498 !$OMP )
3499 
3500 !$OMP SINGLE
3501 !$ ALLOCATE (locks(number_of_locks))
3502 !$OMP END SINGLE
3503 
3504 !$OMP DO
3505 !$ DO lock = 1, number_of_locks
3506 !$ call omp_init_lock(locks(lock))
3507 !$ END DO
3508 !$OMP END DO
3509 
3510  mepos = 0
3511 !$ mepos = omp_get_thread_num()
3512  DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
3513 
3514  CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
3515  IF (ghost(ikind) .OR. ghost(jkind) .OR. floating(ikind) .OR. floating(jkind)) cycle
3516  IF (dispersion_env%nd3_exclude_pair > 0) THEN
3517  CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, exclude=exclude_pair)
3518  IF (exclude_pair) cycle
3519  END IF
3520 
3521  rcc = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
3522  IF (rcc > 1.e-6_dp) THEN
3523  za = atomnumber(ikind)
3524  zb = atomnumber(jkind)
3525  rcova = dispersion_env%rcov(za)
3526  rcovb = dispersion_env%rcov(zb)
3527  CALL cnparam_d3(rcc, rcova, rcovb, dispersion_env%k1, cnab, dcnab)
3528  IF (cnab > eps_cn) THEN
3529 !$OMP ATOMIC
3530  cnumbers(iatom) = cnumbers(iatom) + cnab
3531  IF (iatom /= jatom) THEN
3532 !$OMP ATOMIC
3533  cnumbers(jatom) = cnumbers(jatom) + cnab
3534  END IF
3535  END IF
3536  IF (calculate_forces .OR. debugall .AND. cnab > eps_cn) THEN
3537 !$ CALL omp_set_lock(locks(iatom))
3538  dcnum(iatom)%neighbors = dcnum(iatom)%neighbors + 1
3539  ni = dcnum(iatom)%neighbors
3540  IF (ni > SIZE(dcnum(iatom)%nlist)) THEN
3541  CALL reallocate(dcnum(iatom)%nlist, 1, 2*ni)
3542  CALL reallocate(dcnum(iatom)%dvals, 1, 2*ni)
3543  CALL reallocate(dcnum(iatom)%rik, 1, 3, 1, 2*ni)
3544  END IF
3545  dcnum(iatom)%nlist(ni) = jatom
3546  dcnum(iatom)%dvals(ni) = dcnab
3547  dcnum(iatom)%rik(1:3, ni) = rij(1:3)
3548 !$ CALL omp_unset_lock(locks(iatom))
3549 
3550  IF (iatom /= jatom) THEN
3551 !$ CALL omp_set_lock(locks(jatom))
3552  dcnum(jatom)%neighbors = dcnum(jatom)%neighbors + 1
3553  nj = dcnum(jatom)%neighbors
3554  IF (nj > SIZE(dcnum(jatom)%dvals)) THEN
3555  CALL reallocate(dcnum(jatom)%nlist, 1, 2*nj)
3556  CALL reallocate(dcnum(jatom)%dvals, 1, 2*nj)
3557  CALL reallocate(dcnum(jatom)%rik, 1, 3, 1, 2*nj)
3558  END IF
3559  dcnum(jatom)%nlist(nj) = iatom
3560  dcnum(jatom)%dvals(nj) = dcnab
3561  dcnum(jatom)%rik(1:3, nj) = -rij(1:3)
3562 !$ CALL omp_unset_lock(locks(jatom))
3563  END IF
3564  END IF
3565  END IF
3566  END DO
3567 
3568 !$OMP BARRIER ! Wait for all threads to finish the loop before locks can be freed
3569 !$OMP DO
3570 !$ DO lock = 1, number_of_locks
3571 !$ call omp_destroy_lock(locks(lock))
3572 !$ END DO
3573 !$OMP END DO
3574 !$OMP SINGLE
3575 !$ DEALLOCATE (locks)
3576 !$OMP END SINGLE NOWAIT
3577 
3578 !$OMP END PARALLEL
3579  CALL neighbor_list_iterator_release(nl_iterator)
3580 
3581  CALL timestop(handle)
3582 
3583  END SUBROUTINE d3_cnumber
3584 
3585 ! **************************************************************************************************
3586 
3587 ! **************************************************************************************************
3588 !> \brief ...
3589 !> \param exclude_list List of kind pairs to exclude
3590 !> \param za first kind
3591 !> \param zb second kind
3592 !> \param zc third kind in case of the three-body term
3593 !> \param exclude whether to exclude the interaction or not
3594 ! **************************************************************************************************
3595  SUBROUTINE exclude_d3_kind_pair(exclude_list, za, zb, zc, exclude)
3596 
3597  INTEGER, DIMENSION(:, :), INTENT(IN) :: exclude_list
3598  INTEGER, INTENT(in) :: za, zb
3599  INTEGER, INTENT(in), OPTIONAL :: zc
3600  LOGICAL, INTENT(out) :: exclude
3601 
3602  CHARACTER(LEN=*), PARAMETER :: routinen = 'exclude_d3_kind_pair'
3603 
3604  INTEGER :: handle, i
3605 
3606  CALL timeset(routinen, handle)
3607  exclude = .false.
3608  DO i = 1, SIZE(exclude_list, 1)
3609  IF (exclude_list(i, 1) == za .AND. exclude_list(i, 2) == zb .OR. &
3610  exclude_list(i, 1) == zb .AND. exclude_list(i, 2) == za) THEN
3611  exclude = .true.
3612  EXIT
3613  END IF
3614  IF (PRESENT(zc)) THEN
3615  IF (exclude_list(i, 1) == za .AND. exclude_list(i, 2) == zc .OR. &
3616  exclude_list(i, 1) == zc .AND. exclude_list(i, 2) == za .OR. &
3617  exclude_list(i, 1) == zb .AND. exclude_list(i, 2) == zc .OR. &
3618  exclude_list(i, 1) == zc .AND. exclude_list(i, 2) == zb) THEN
3619  exclude = .true.
3620  EXIT
3621  END IF
3622  END IF
3623  END DO
3624 
3625  CALL timestop(handle)
3626 
3627  END SUBROUTINE exclude_d3_kind_pair
3628 
3629 END MODULE qs_dispersion_pairpot
hash
static unsigned int hash(const dbm_task_t task)
Private hash function based on Szudzik's elegant pairing. Using unsigned int to return a positive num...
Definition: dbm_multiply_cpu.c:57
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
uppercase
subroutine uppercase(string)
...
Definition: dumpdcd.F:1376
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition: atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition: atprop_types.F:98
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::grimme2006
integer, save, public grimme2006
Definition: bibliography.F:43
bibliography::goerigk2017
integer, save, public goerigk2017
Definition: bibliography.F:43
bibliography::grimme2011
integer, save, public grimme2011
Definition: bibliography.F:43
bibliography::grimme2010
integer, save, public grimme2010
Definition: bibliography.F:43
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
cell_types::plane_distance
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition: cell_types.F:252
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Definition: cp_parser_methods.F:421
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition: cp_parser_types.F:109
cp_parser_types::parser_create
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Definition: cp_parser_types.F:147
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::xc_vdw_fun_none
integer, parameter, public xc_vdw_fun_none
Definition: input_constants.F:592
input_constants::vdw_pairpot_dftd3
integer, parameter, public vdw_pairpot_dftd3
Definition: input_constants.F:595
input_constants::vdw_pairpot_dftd2
integer, parameter, public vdw_pairpot_dftd2
Definition: input_constants.F:595
input_constants::xc_vdw_fun_pairpot
integer, parameter, public xc_vdw_fun_pairpot
Definition: input_constants.F:592
input_constants::vdw_pairpot_dftd3bj
integer, parameter, public vdw_pairpot_dftd3bj
Definition: input_constants.F:595
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::kcalmol
real(kind=dp), parameter, public kcalmol
Definition: physcon.F:171
physcon::kjmol
real(kind=dp), parameter, public kjmol
Definition: physcon.F:168
physcon::bohr
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
qs_dispersion_pairpot
Calculation of dispersion using pair potentials.
Definition: qs_dispersion_pairpot.F:12
qs_dispersion_pairpot::qs_scaling_dftd3bj
subroutine, public qs_scaling_dftd3bj(s6, a1, s8, a2, vdw_section)
...
Definition: qs_dispersion_pairpot.F:828
qs_dispersion_pairpot::qs_dispersion_pairpot_init
subroutine, public qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
...
Definition: qs_dispersion_pairpot.F:116
qs_dispersion_pairpot::dcnum_distribute
subroutine, public dcnum_distribute(dcnum, para_env)
...
Definition: qs_dispersion_pairpot.F:3369
qs_dispersion_pairpot::d3_cnumber
subroutine, public d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, calculate_forces, debugall)
...
Definition: qs_dispersion_pairpot.F:3449
qs_dispersion_pairpot::qs_dispersion_setcn
subroutine, public qs_dispersion_setcn(atomic_kind_set, qs_kind_set, dispersion_env, para_env)
...
Definition: qs_dispersion_pairpot.F:359
qs_dispersion_pairpot::qs_scaling_init
subroutine, public qs_scaling_init(scaling, vdw_section)
...
Definition: qs_dispersion_pairpot.F:413
qs_dispersion_pairpot::calculate_dispersion_pairpot
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition: qs_dispersion_pairpot.F:1245
qs_dispersion_pairpot::qs_scaling_dftd3
subroutine, public qs_scaling_dftd3(s6, sr6, s8, vdw_section)
...
Definition: qs_dispersion_pairpot.F:454
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition: qs_dispersion_types.F:12
qs_dispersion_types::dftd2_pp
integer, parameter, public dftd2_pp
Definition: qs_dispersion_types.F:28
qs_dispersion_types::dftd3_pp
integer, parameter, public dftd3_pp
Definition: qs_dispersion_types.F:29
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::set_qs_kind
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Definition: qs_kind_types.F:2799
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition: string_utilities.F:3361
virial_methods
Definition: virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition: virial_methods.F:139
virial_types
Definition: virial_types.F:13