(git:6344619)
Loading...
Searching...
No Matches
qs_dftb_parameters.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \author JGH (27.02.2007)
10! **************************************************************************************************
11MODULE qs_dftb_parameters
12
13 USE atomic_kind_types, ONLY: atomic_kind_type,&
14 get_atomic_kind
15 USE cp_control_types, ONLY: dftb_control_type
16 USE cp_files, ONLY: close_file,&
17 open_file
18 USE cp_log_handling, ONLY: cp_get_default_logger,&
19 cp_logger_type
20 USE cp_output_handling, ONLY: cp_p_file,&
21 cp_print_key_finished_output,&
22 cp_print_key_should_output,&
23 cp_print_key_unit_nr
24 USE cp_parser_methods, ONLY: parser_get_next_line,&
25 parser_get_object
26 USE cp_parser_types, ONLY: cp_parser_type,&
27 parser_create,&
28 parser_release
29 USE external_potential_types, ONLY: set_potential
30 USE input_constants, ONLY: dispersion_uff
31 USE input_section_types, ONLY: section_vals_type
32 USE kinds, ONLY: default_path_length,&
33 default_string_length,&
34 dp
35 USE mathconstants, ONLY: pi
36 USE message_passing, ONLY: mp_para_env_type
37 USE physcon, ONLY: angstrom,&
38 kcalmol
39 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
40 USE qs_dftb_types, ONLY: qs_dftb_atom_type,&
41 qs_dftb_pairpot_create,&
42 qs_dftb_pairpot_init,&
43 qs_dftb_pairpot_type
44 USE qs_dftb_utils, ONLY: allocate_dftb_atom_param,&
45 get_dftb_atom_param,&
46 set_dftb_atom_param
47 USE qs_kind_types, ONLY: get_qs_kind,&
48 qs_kind_type,&
49 set_qs_kind
50 USE string_utilities, ONLY: uppercase
51#include "../base/base_uses.f90"
52
53 IMPLICIT NONE
54
55 PRIVATE
56
57! *** Global parameters ***
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_parameters'
60
61 REAL(dp), PARAMETER :: slako_d0 = 1._dp
62
63! *** Public subroutines ***
64
65 PUBLIC :: qs_dftb_param_init
66
67CONTAINS
68
69! **************************************************************************************************
70!> \brief ...
71!> \param atomic_kind_set ...
72!> \param qs_kind_set ...
73!> \param dftb_control ...
74!> \param dftb_potential ...
75!> \param subsys_section ...
76!> \param para_env ...
77! **************************************************************************************************
78 SUBROUTINE qs_dftb_param_init(atomic_kind_set, qs_kind_set, dftb_control, dftb_potential, &
79 subsys_section, para_env)
80 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
81 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
82 TYPE(dftb_control_type), INTENT(inout) :: dftb_control
83 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
84 POINTER :: dftb_potential
85 TYPE(section_vals_type), POINTER :: subsys_section
86 TYPE(mp_para_env_type), POINTER :: para_env
87
88 CHARACTER(LEN=2) :: iel, jel
89 CHARACTER(LEN=6) :: cspline
90 CHARACTER(LEN=default_path_length) :: file_name
91 CHARACTER(LEN=default_path_length), ALLOCATABLE, &
92 DIMENSION(:, :) :: sk_files
93 CHARACTER(LEN=default_string_length) :: iname, jname, name_a, name_b, skfn
94 INTEGER :: ikind, isp, jkind, k, l, l1, l2, llm, &
95 lmax, lmax_a, lmax_b, lp, m, n_urpoly, &
96 ngrd, nkind, output_unit, runit, &
97 spdim, z
98 LOGICAL :: at_end, found, ldum, search, sklist
99 REAL(dp) :: da, db, dgrd, dij, energy, eps_disp, ra, &
100 radmax, rb, rcdisp, rmax6, s_cut, xij, &
101 zeff
102 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fmat, scoeff, smat, spxr
103 REAL(dp), DIMENSION(0:3) :: eta, occupation, skself
104 REAL(dp), DIMENSION(10) :: fwork, swork, uwork
105 REAL(dp), DIMENSION(1:2) :: surr
106 REAL(dp), DIMENSION(1:3) :: srep
107 TYPE(cp_logger_type), POINTER :: logger
108 TYPE(qs_dftb_atom_type), POINTER :: dftb_atom_a, dftb_atom_b
109
110 output_unit = -1
111 NULLIFY (logger)
112 logger => cp_get_default_logger()
113 IF (btest(cp_print_key_should_output(logger%iter_info, subsys_section, &
114 "PRINT%KINDS/BASIS_SET"), cp_p_file)) THEN
115 output_unit = cp_print_key_unit_nr(logger, subsys_section, &
116 "PRINT%KINDS", extension=".Log")
117 IF (output_unit > 0) THEN
118 WRITE (output_unit, "(/,A)") " DFTB| A set of relativistic DFTB "// &
119 "parameters for material sciences."
120 WRITE (output_unit, "(A)") " DFTB| J. Frenzel, N. Jardillier, A.F. Oliveira,"// &
121 " T. Heine, G. Seifert"
122 WRITE (output_unit, "(A)") " DFTB| TU Dresden, 2002-2007"
123 WRITE (output_unit, "(/,A)") " DFTB| Non-SCC parameters "
124 WRITE (output_unit, "(A,T25,A)") " DFTB| C,H :", &
125 " D. Porezag et al, PRB 51 12947 (1995)"
126 WRITE (output_unit, "(A,T25,A)") " DFTB| B,N :", &
127 " J. Widany et al, PRB 53 4443 (1996)"
128 WRITE (output_unit, "(A,T25,A)") " DFTB| Li,Na,K,Cl :", &
129 " S. Hazebroucq et al, JCP 123 134510 (2005)"
130 WRITE (output_unit, "(A,T25,A)") " DFTB| F :", &
131 " T. Heine et al, JCSoc-Perkins Trans 2 707 (1999)"
132 WRITE (output_unit, "(A,T25,A)") " DFTB| Mo,S :", &
133 " G. Seifert et al, PRL 85 146 (2000)"
134 WRITE (output_unit, "(A,T25,A)") " DFTB| P :", &
135 " G. Seifert et al, EPS 16 341 (2001)"
136 WRITE (output_unit, "(A,T25,A)") " DFTB| Sc,N,C :", &
137 " M. Krause et al, JCP 115 6596 (2001)"
138 END IF
139 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
140 "PRINT%KINDS")
141 END IF
142
143 sklist = (dftb_control%sk_file_list /= "")
144
145 nkind = SIZE(atomic_kind_set)
146 ALLOCATE (sk_files(nkind, nkind))
147 ! allocate potential structures
148 ALLOCATE (dftb_potential(nkind, nkind))
149 CALL qs_dftb_pairpot_init(dftb_potential)
150
151 DO ikind = 1, nkind
152 CALL get_atomic_kind(atomic_kind_set(ikind), name=iname, element_symbol=iel)
153 CALL uppercase(iname)
154 CALL uppercase(iel)
155 ldum = qmmm_ff_precond_only_qm(iname)
156 DO jkind = 1, nkind
157 CALL get_atomic_kind(atomic_kind_set(jkind), name=jname, element_symbol=jel)
158 CALL uppercase(jname)
159 CALL uppercase(jel)
160 ldum = qmmm_ff_precond_only_qm(jname)
161 found = .false.
162 DO k = 1, SIZE(dftb_control%sk_pair_list, 2)
163 name_a = trim(dftb_control%sk_pair_list(1, k))
164 name_b = trim(dftb_control%sk_pair_list(2, k))
165 CALL uppercase(name_a)
166 CALL uppercase(name_b)
167 IF ((iname == name_a .AND. jname == name_b)) THEN
168 sk_files(ikind, jkind) = trim(dftb_control%sk_file_path)//"/"// &
169 trim(dftb_control%sk_pair_list(3, k))
170 found = .true.
171 EXIT
172 END IF
173 END DO
174 IF (.NOT. found .AND. sklist) THEN
175 file_name = trim(dftb_control%sk_file_path)//"/"// &
176 trim(dftb_control%sk_file_list)
177 block
178 TYPE(cp_parser_type) :: parser
179 CALL parser_create(parser, file_name, para_env=para_env)
180 DO
181 at_end = .false.
182 CALL parser_get_next_line(parser, 1, at_end)
183 IF (at_end) EXIT
184 CALL parser_get_object(parser, name_a, lower_to_upper=.true.)
185 CALL parser_get_object(parser, name_b, lower_to_upper=.true.)
186 !Checking Names
187 IF ((iname == name_a .AND. jname == name_b)) THEN
188 CALL parser_get_object(parser, skfn, string_length=8)
189 sk_files(ikind, jkind) = trim(dftb_control%sk_file_path)//"/"// &
190 trim(skfn)
191 found = .true.
192 EXIT
193 END IF
194 !Checking Element
195 IF ((iel == name_a .AND. jel == name_b)) THEN
196 CALL parser_get_object(parser, skfn, string_length=8)
197 sk_files(ikind, jkind) = trim(dftb_control%sk_file_path)//"/"// &
198 trim(skfn)
199 found = .true.
200 EXIT
201 END IF
202 END DO
203 CALL parser_release(parser)
204 END block
205 END IF
206 IF (.NOT. found) THEN
207 CALL cp_abort(__location__, &
208 "Failure in assigning KINDS <"//trim(iname)//"> and <"//trim(jname)// &
209 "> to a DFTB interaction pair!")
210 END IF
211 END DO
212 END DO
213 ! reading the files
214 ! read all pairs, equal kind first
215 DO ikind = 1, nkind
216 CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
217
218 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
219 IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
220 CALL allocate_dftb_atom_param(dftb_atom_a)
221 CALL set_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
222 END IF
223
224 ! read all pairs, equal kind first
225 jkind = ikind
226
227 CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
228 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
229
230 IF (output_unit > 0) THEN
231 WRITE (output_unit, "(A,T30,A50)") " DFTB| Reading parameter file ", &
232 adjustr(trim(sk_files(jkind, ikind)))
233 END IF
234 skself = 0._dp
235 eta = 0._dp
236 occupation = 0._dp
237 IF (para_env%is_source()) THEN
238 CALL open_file(file_name=sk_files(jkind, ikind), unit_number=runit)
239 ! grid density and number of grid poin ts
240 READ (runit, fmt=*, END=1, err=1) dgrd, ngrd
241!
242! ngrd -1 ?
243! In Slako tables, the grid starts at 0.0, in deMon it starts with dgrd
244!
245 ngrd = ngrd - 1
246!
247 ! orbital energy, total energy, hardness, occupation
248 READ (runit, fmt=*, END=1, err=1) skself(2:0:-1), energy, &
249 eta(2:0:-1), occupation(2:0:-1)
250 ! repulsive potential as polynomial
251 READ (runit, fmt=*, END=1, err=1) uwork(1:10)
252 n_urpoly = 0
253 IF (dot_product(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
254 n_urpoly = 1
255 DO k = 2, 9
256 IF (abs(uwork(k)) >= 1.e-12_dp) n_urpoly = k
257 END DO
258 END IF
259! Polynomials of length 1 are not allowed, it seems we should use spline after all
260! This is creative guessing!
261 IF (n_urpoly < 2) n_urpoly = 0
262 END IF
263
264 CALL para_env%bcast(n_urpoly)
265 CALL para_env%bcast(uwork)
266 CALL para_env%bcast(ngrd)
267 CALL para_env%bcast(dgrd)
268
269 CALL para_env%bcast(skself)
270 CALL para_env%bcast(energy)
271 CALL para_env%bcast(eta)
272 CALL para_env%bcast(occupation)
273
274 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, &
275 z=z, zeff=sum(occupation), defined=.true., &
276 skself=skself, energy=energy, eta=eta, occupation=occupation)
277
278 ! Slater-Koster table
279 ALLOCATE (fmat(ngrd, 10))
280 ALLOCATE (smat(ngrd, 10))
281 IF (para_env%is_source()) THEN
282 DO k = 1, ngrd
283 READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
284 fmat(k, 1:10) = fwork(1:10)
285 smat(k, 1:10) = swork(1:10)
286 END DO
287 END IF
288 CALL para_env%bcast(fmat)
289 CALL para_env%bcast(smat)
290
291 !
292 ! Determine lmax for atom type.
293 ! An atomic orbital is 'active' if either its onsite energy is different from zero,
294 ! or
295 ! if this matrix element contains non-zero elements.
296 ! The sigma interactions are sufficient for that.
297 ! In the DFTB-Slako convention they are on orbital 10 (s-s-sigma),
298 ! 7 (p-p-sigma) and 3 (d-d-sigma).
299 !
300 ! We also allow lmax to be set in the input (in KIND)
301 !
302 CALL get_qs_kind(qs_kind_set(ikind), lmax_dftb=lmax)
303 IF (lmax < 0) THEN
304 lmax = 0
305 DO l = 0, 3
306 SELECT CASE (l)
307 CASE DEFAULT
308 cpabort("Only 0, 1, 2 are supported as the value of l")
309 CASE (0)
310 lp = 10
311 CASE (1)
312 lp = 7
313 CASE (2)
314 lp = 3
315 CASE (3)
316 lp = 3 ! this is wrong but we don't allow f anyway
317 END SELECT
318 ! Technical note: In some slako files dummies are included in the
319 ! first matrix elements, so remove them.
320 IF ((abs(skself(l)) > 0._dp) .OR. &
321 (sum(abs(fmat(ngrd/10:ngrd, lp))) > 0._dp)) lmax = l
322 END DO
323 ! l=2 (d) is maximum
324 lmax = min(2, lmax)
325 END IF
326 IF (lmax > 2) THEN
327 CALL cp_abort(__location__, "Maximum L allowed is d. "// &
328 "Use KIND/LMAX_DFTB to set smaller values if needed.")
329 END IF
330 !
331 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, &
332 lmax=lmax, natorb=(lmax + 1)**2)
333
334 spdim = 0
335 IF (n_urpoly == 0) THEN
336 IF (para_env%is_source()) THEN
337 ! Look for spline representation of repulsive potential
338 search = .true.
339 DO WHILE (search)
340 READ (runit, fmt='(A6)', END=1, err=1) cspline
341 IF (cspline == 'Spline') THEN
342 search = .false.
343 ! spline dimension and left-hand cutoff
344 READ (runit, fmt=*, END=1, err=1) spdim, s_cut
345 ALLOCATE (spxr(spdim, 2))
346 ALLOCATE (scoeff(spdim, 4))
347 ! e-functions describing left-hand extrapolation
348 READ (runit, fmt=*, END=1, err=1) srep(1:3)
349 DO isp = 1, spdim - 1
350 ! location and coefficients of 'normal' spline range
351 READ (runit, fmt=*, END=1, err=1) spxr(isp, 1:2), scoeff(isp, 1:4)
352 END DO
353 ! last point has 2 more coefficients
354 READ (runit, fmt=*, END=1, err=1) spxr(spdim, 1:2), scoeff(spdim, 1:4), surr(1:2)
355 END IF
356 END DO
357 END IF
358 END IF
359
360 IF (para_env%is_source()) THEN
361 CALL close_file(unit_number=runit)
362 END IF
363
364 CALL para_env%bcast(spdim)
365 IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
366 ALLOCATE (spxr(spdim, 2))
367 ALLOCATE (scoeff(spdim, 4))
368 END IF
369 IF (spdim > 0) THEN
370 CALL para_env%bcast(spxr)
371 CALL para_env%bcast(scoeff)
372 CALL para_env%bcast(surr)
373 CALL para_env%bcast(srep)
374 CALL para_env%bcast(s_cut)
375 END IF
376
377 ! store potential data
378 ! allocate data
379 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
380 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
381 llm = 0
382 DO l1 = 0, max(lmax_a, lmax_b)
383 DO l2 = 0, min(l1, lmax_a, lmax_b)
384 DO m = 0, l2
385 llm = llm + 1
386 END DO
387 END DO
388 END DO
389 CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
390 ngrd, llm, spdim)
391
392 ! repulsive potential
393 dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
394 dftb_potential(ikind, jkind)%urep_cut = uwork(10)
395 dftb_potential(ikind, jkind)%urep(:) = 0._dp
396 dftb_potential(ikind, jkind)%urep(1) = uwork(10)
397 dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
398
399 ! Slater-Koster tables
400 dftb_potential(ikind, jkind)%dgrd = dgrd
401 CALL skreorder(fmat, lmax_a, lmax_b)
402 dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
403 CALL skreorder(smat, lmax_a, lmax_b)
404 dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
405 dftb_potential(ikind, jkind)%ngrdcut = ngrd + int(slako_d0/dgrd)
406 ! Splines
407 IF (spdim > 0) THEN
408 dftb_potential(ikind, jkind)%s_cut = s_cut
409 dftb_potential(ikind, jkind)%srep = srep
410 dftb_potential(ikind, jkind)%spxr = spxr
411 dftb_potential(ikind, jkind)%scoeff = scoeff
412 dftb_potential(ikind, jkind)%surr = surr
413 END IF
414
415 DEALLOCATE (fmat)
416 DEALLOCATE (smat)
417 IF (spdim > 0) THEN
418 DEALLOCATE (spxr)
419 DEALLOCATE (scoeff)
420 END IF
421
422 END DO
423
424 ! no all other pairs
425 DO ikind = 1, nkind
426 CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
427 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
428
429 IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
430 CALL allocate_dftb_atom_param(dftb_atom_a)
431 CALL set_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
432 END IF
433
434 DO jkind = 1, nkind
435
436 IF (ikind == jkind) cycle
437 CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
438 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
439
440 IF (output_unit > 0) THEN
441 WRITE (output_unit, "(A,T30,A50)") " DFTB| Reading parameter file ", &
442 adjustr(trim(sk_files(ikind, jkind)))
443 END IF
444 skself = 0._dp
445 eta = 0._dp
446 occupation = 0._dp
447 IF (para_env%is_source()) THEN
448 CALL open_file(file_name=sk_files(ikind, jkind), unit_number=runit)
449 ! grid density and number of grid poin ts
450 READ (runit, fmt=*, END=1, err=1) dgrd, ngrd
451!
452! ngrd -1 ?
453! In Slako tables, the grid starts at 0.0, in deMon it starts with dgrd
454!
455 ngrd = ngrd - 1
456!
457 IF (ikind == jkind) THEN
458 ! orbital energy, total energy, hardness, occupation
459 READ (runit, fmt=*, END=1, err=1) skself(2:0:-1), energy, &
460 eta(2:0:-1), occupation(2:0:-1)
461 END IF
462 ! repulsive potential as polynomial
463 READ (runit, fmt=*, END=1, err=1) uwork(1:10)
464 n_urpoly = 0
465 IF (dot_product(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
466 n_urpoly = 1
467 DO k = 2, 9
468 IF (abs(uwork(k)) >= 1.e-12_dp) n_urpoly = k
469 END DO
470 END IF
471! Polynomials of length 1 are not allowed, it seems we should use spline after all
472! This is creative guessing!
473 IF (n_urpoly < 2) n_urpoly = 0
474 END IF
475
476 CALL para_env%bcast(n_urpoly)
477 CALL para_env%bcast(uwork)
478 CALL para_env%bcast(ngrd)
479 CALL para_env%bcast(dgrd)
480
481 ! Slater-Koster table
482 ALLOCATE (fmat(ngrd, 10))
483 ALLOCATE (smat(ngrd, 10))
484 IF (para_env%is_source()) THEN
485 DO k = 1, ngrd
486 READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
487 fmat(k, 1:10) = fwork(1:10)
488 smat(k, 1:10) = swork(1:10)
489 END DO
490 END IF
491 CALL para_env%bcast(fmat)
492 CALL para_env%bcast(smat)
493
494 spdim = 0
495 IF (n_urpoly == 0) THEN
496 IF (para_env%is_source()) THEN
497 ! Look for spline representation of repulsive potential
498 search = .true.
499 DO WHILE (search)
500 READ (runit, fmt='(A6)', END=1, err=1) cspline
501 IF (cspline == 'Spline') THEN
502 search = .false.
503 ! spline dimension and left-hand cutoff
504 READ (runit, fmt=*, END=1, err=1) spdim, s_cut
505 ALLOCATE (spxr(spdim, 2))
506 ALLOCATE (scoeff(spdim, 4))
507 ! e-functions describing left-hand extrapolation
508 READ (runit, fmt=*, END=1, err=1) srep(1:3)
509 DO isp = 1, spdim - 1
510 ! location and coefficients of 'normal' spline range
511 READ (runit, fmt=*, END=1, err=1) spxr(isp, 1:2), scoeff(isp, 1:4)
512 END DO
513 ! last point has 2 more coefficients
514 READ (runit, fmt=*, END=1, err=1) spxr(spdim, 1:2), scoeff(spdim, 1:4), surr(1:2)
515 END IF
516 END DO
517 END IF
518 END IF
519
520 IF (para_env%is_source()) THEN
521 CALL close_file(unit_number=runit)
522 END IF
523
524 CALL para_env%bcast(spdim)
525 IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
526 ALLOCATE (spxr(spdim, 2))
527 ALLOCATE (scoeff(spdim, 4))
528 END IF
529 IF (spdim > 0) THEN
530 CALL para_env%bcast(spxr)
531 CALL para_env%bcast(scoeff)
532 CALL para_env%bcast(surr)
533 CALL para_env%bcast(srep)
534 CALL para_env%bcast(s_cut)
535 END IF
536
537 ! store potential data
538 ! allocate data
539 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
540 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
541 llm = 0
542 DO l1 = 0, max(lmax_a, lmax_b)
543 DO l2 = 0, min(l1, lmax_a, lmax_b)
544 DO m = 0, l2
545 llm = llm + 1
546 END DO
547 END DO
548 END DO
549 CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
550 ngrd, llm, spdim)
551
552 ! repulsive potential
553 dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
554 dftb_potential(ikind, jkind)%urep_cut = uwork(10)
555 dftb_potential(ikind, jkind)%urep(:) = 0._dp
556 dftb_potential(ikind, jkind)%urep(1) = uwork(10)
557 dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
558
559 ! Slater-Koster tables
560 dftb_potential(ikind, jkind)%dgrd = dgrd
561 CALL skreorder(fmat, lmax_a, lmax_b)
562 dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
563 CALL skreorder(smat, lmax_a, lmax_b)
564 dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
565 dftb_potential(ikind, jkind)%ngrdcut = ngrd + int(slako_d0/dgrd)
566 ! Splines
567 IF (spdim > 0) THEN
568 dftb_potential(ikind, jkind)%s_cut = s_cut
569 dftb_potential(ikind, jkind)%srep = srep
570 dftb_potential(ikind, jkind)%spxr = spxr
571 dftb_potential(ikind, jkind)%scoeff = scoeff
572 dftb_potential(ikind, jkind)%surr = surr
573 END IF
574
575 DEALLOCATE (fmat)
576 DEALLOCATE (smat)
577 IF (spdim > 0) THEN
578 DEALLOCATE (spxr)
579 DEALLOCATE (scoeff)
580 END IF
581
582 END DO
583 END DO
584
585 DEALLOCATE (sk_files)
586
587 ! read dispersion parameters (UFF type)
588 IF (dftb_control%dispersion) THEN
589
590 IF (dftb_control%dispersion_type == dispersion_uff) THEN
591 file_name = trim(dftb_control%sk_file_path)//"/"// &
592 trim(dftb_control%uff_force_field)
593 block
594 TYPE(cp_parser_type) :: parser
595 DO ikind = 1, nkind
596 CALL get_atomic_kind(atomic_kind_set(ikind), name=iname)
597 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
598
599 m = len_trim(iname)
600 CALL parser_create(parser, file_name, para_env=para_env)
601 found = .false.
602 DO
603 at_end = .false.
604 CALL parser_get_next_line(parser, 1, at_end)
605 IF (at_end) EXIT
606 CALL parser_get_object(parser, name_a)
607 ! parser is no longer removing leading quotes
608 IF (name_a(1:1) == '"') name_a(1:m) = name_a(2:m + 1)
609 IF (name_a(1:m) == trim(iname)) THEN
610 CALL parser_get_object(parser, rb)
611 CALL parser_get_object(parser, rb)
612 CALL parser_get_object(parser, ra)
613 CALL parser_get_object(parser, da)
614 found = .true.
615 ra = ra/angstrom
616 da = da/kcalmol
617 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, name=iname, xi=ra, di=da)
618 EXIT
619 END IF
620 END DO
621 CALL parser_release(parser)
622 END DO
623 END block
624 END IF
625
626 END IF
627
628 ! extract simple atom interaction radii
629 DO ikind = 1, nkind
630 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
631 radmax = (dftb_potential(ikind, ikind)%ngrdcut + 1)* &
632 dftb_potential(ikind, ikind)%dgrd*0.5_dp
633 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=radmax)
634 END DO
635 DO ikind = 1, nkind
636 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
637 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
638 DO jkind = 1, nkind
639 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
640 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, cutoff=rb)
641 radmax = (dftb_potential(ikind, jkind)%ngrdcut + 1)* &
642 dftb_potential(ikind, jkind)%dgrd
643 IF (ra + rb < radmax) THEN
644 ra = ra + (radmax - ra - rb)*0.5_dp
645 rb = rb + (radmax - ra - rb)*0.5_dp
646 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
647 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_b, cutoff=rb)
648 END IF
649 END DO
650 END DO
651
652 ! set correct core charge in potential
653 DO ikind = 1, nkind
654 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
655 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, zeff=zeff)
656 CALL set_potential(potential=qs_kind_set(ikind)%all_potential, &
657 zeff=zeff, zeff_correction=0.0_dp)
658 END DO
659
660 ! setup DFTB3 parameters
661 IF (dftb_control%dftb3_diagonal) THEN
662 DO ikind = 1, nkind
663 CALL get_qs_kind(qs_kind_set(ikind), dftb3_param=db)
664 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
665 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, dudq=db)
666 END DO
667 END IF
668
669 ! setup dispersion parameters (UFF type)
670 IF (dftb_control%dispersion) THEN
671 IF (dftb_control%dispersion_type == dispersion_uff) THEN
672 eps_disp = dftb_control%eps_disp
673 DO ikind = 1, nkind
674 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
675 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, xi=ra, di=da)
676 rcdisp = 0._dp
677 DO jkind = 1, nkind
678 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
679 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, xi=rb, di=db)
680 xij = sqrt(ra*rb)
681 dij = sqrt(da*db)
682 dftb_potential(ikind, jkind)%xij = xij
683 dftb_potential(ikind, jkind)%dij = dij
684 dftb_potential(ikind, jkind)%x0ij = xij*(0.5_dp**(1.0_dp/6.0_dp))
685 dftb_potential(ikind, jkind)%a = dij*396.0_dp/25.0_dp
686 dftb_potential(ikind, jkind)%b = &
687 dij/(xij**5)*672.0_dp*2.0_dp**(5.0_dp/6.0_dp)/25.0_dp
688 dftb_potential(ikind, jkind)%c = &
689 -dij/(xij**10)*2.0_dp**(2.0_dp/3.0_dp)*552.0_dp/25.0_dp
690 rmax6 = ((8._dp*pi*dij/eps_disp)*xij**6)**0.25_dp
691 rcdisp = max(rcdisp, rmax6*0.5_dp)
692 END DO
693 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, rcdisp=rcdisp)
694 END DO
695 END IF
696 END IF
697
698 RETURN
699
7001 CONTINUE
701 ! Many instances of READ (..., END=1, err=1) gets conflated here
702 ! TODO: overhaul the file parser to handle errors separately
703 cpabort("Something went wrong while reading DFTB parameter file")
704
705 END SUBROUTINE qs_dftb_param_init
706
707! **************************************************************************************************
708!> \brief Transform Slako format in l1/l2/m format
709!> \param xmat ...
710!> \param la ...
711!> \param lb ...
712!> \par Notes
713!> Slako tables from Dresden/Paderborn/Heidelberg groups are
714!> stored in the following native format:
715!>
716!> Convention: Higher angular momenta are always on the right-hand side
717!>
718!> 1: d - d - delta
719!> 2: d - d - pi
720!> 3: d - d - sigma
721!> 4: p - d - pi
722!> 5: p - d - sigma
723!> 6: p - p - pi
724!> 7: p - p - sigma
725!> 8: d - s - sigma
726!> 9: p - s - sigma
727!> 10: s - s - sigma
728!> \version 1.0
729! **************************************************************************************************
730 SUBROUTINE skreorder(xmat, la, lb)
731 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: xmat
732 INTEGER, INTENT(IN) :: la, lb
733
734 INTEGER :: i, l1, l2, llm, m
735 REAL(dp) :: skllm(0:3, 0:3, 0:3)
736
737 DO i = 1, SIZE(xmat, 1)
738 skllm = 0._dp
739 skllm(0, 0, 0) = xmat(i, 10)
740 skllm(1, 0, 0) = xmat(i, 9)
741 skllm(2, 0, 0) = xmat(i, 8)
742 skllm(1, 1, 1) = xmat(i, 7)
743 skllm(1, 1, 0) = xmat(i, 6)
744 skllm(2, 1, 1) = xmat(i, 5)
745 skllm(2, 1, 0) = xmat(i, 4)
746 skllm(2, 2, 2) = xmat(i, 3)
747 skllm(2, 2, 1) = xmat(i, 2)
748 skllm(2, 2, 0) = xmat(i, 1)
749 llm = 0
750 DO l1 = 0, max(la, lb)
751 DO l2 = 0, min(l1, la, lb)
752 DO m = 0, l2
753 llm = llm + 1
754 xmat(i, llm) = skllm(l1, l2, m)
755 END DO
756 END DO
757 END DO
758 END DO
759 !
760 END SUBROUTINE skreorder
761
762END MODULE qs_dftb_parameters
763
cp_parser_methods::parser_get_object
Definition cp_parser_methods.F:54
external_potential_types::set_potential
Definition external_potential_types.F:310
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:838
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1058
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:104
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:330
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Definition cp_parser_methods.F:421
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition cp_parser_types.F:109
cp_parser_types::parser_create
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Definition cp_parser_types.F:147
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::dispersion_uff
integer, parameter, public dispersion_uff
Definition input_constants.F:1286
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::kcalmol
real(kind=dp), parameter, public kcalmol
Definition physcon.F:171
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
qmmm_ff_fist
Definition qmmm_ff_fist.F:11
qmmm_ff_fist::qmmm_ff_precond_only_qm
logical function, public qmmm_ff_precond_only_qm(id1, id2, id3, id4, is_link)
This function handles the atom names and modifies the "_QM_" prefix, in order to find the parameters ...
Definition qmmm_ff_fist.F:39
qs_dftb_parameters
Definition qs_dftb_parameters.F:11
qs_dftb_parameters::qs_dftb_param_init
subroutine, public qs_dftb_param_init(atomic_kind_set, qs_kind_set, dftb_control, dftb_potential, subsys_section, para_env)
...
Definition qs_dftb_parameters.F:80
qs_dftb_types
Definition of the DFTB parameter types.
Definition qs_dftb_types.F:12
qs_dftb_types::qs_dftb_pairpot_init
subroutine, public qs_dftb_pairpot_init(pairpot)
...
Definition qs_dftb_types.F:83
qs_dftb_types::qs_dftb_pairpot_create
subroutine, public qs_dftb_pairpot_create(pairpot, ngrd, llm, spdim)
...
Definition qs_dftb_types.F:107
qs_dftb_utils
Working with the DFTB parameter types.
Definition qs_dftb_utils.F:12
qs_dftb_utils::set_dftb_atom_param
subroutine, public set_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition qs_dftb_utils.F:178
qs_dftb_utils::get_dftb_atom_param
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition qs_dftb_utils.F:124
qs_dftb_utils::allocate_dftb_atom_param
subroutine, public allocate_dftb_atom_param(dftb_parameter)
...
Definition qs_dftb_utils.F:62
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_kind_types::set_qs_kind
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, hund_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Definition qs_kind_types.F:3328
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_control_types::dftb_control_type
Definition cp_control_types.F:152
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_parser_types::cp_parser_type
represent a parser
Definition cp_parser_types.F:77
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
qs_dftb_types::qs_dftb_atom_type
Definition qs_dftb_types.F:27
qs_dftb_types::qs_dftb_pairpot_type
Definition qs_dftb_types.F:47
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190