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