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