(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
12MODULE 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
105CONTAINS
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
3629END 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
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:48
cell_types::pbc
Definition cell_types.F:92
cp_parser_methods::parser_get_object
Definition cp_parser_methods.F:54
memory_utilities::reallocate
Definition memory_utilities.F:27
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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_parser_types::cp_parser_type
represent a parser
Definition cp_parser_types.F:77
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_pairpot::dcnum_type
Definition qs_dispersion_pairpot.F:94
qs_dispersion_types::qs_atom_dispersion_type
Definition qs_dispersion_types.F:99
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:31
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26