(git:d18deda)
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-2025 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) &
208 CALL cp_abort(__location__, &
209 "Failure in assigning KINDS <"//trim(iname)//"> and <"//trim(jname)// &
210 "> to a DFTB interaction pair!")
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 runit = get_unit_number()
239 CALL open_file(file_name=sk_files(jkind, ikind), unit_number=runit)
240 ! grid density and number of grid poin ts
241 READ (runit, fmt=*, END=1, err=1) dgrd, ngrd
242!
243! ngrd -1 ?
244! In Slako tables, the grid starts at 0.0, in deMon it starts with dgrd
245!
246 ngrd = ngrd - 1
247!
248 ! orbital energy, total energy, hardness, occupation
249 READ (runit, fmt=*, END=1, err=1) skself(2:0:-1), energy, &
250 eta(2:0:-1), occupation(2:0:-1)
251 ! repulsive potential as polynomial
252 READ (runit, fmt=*, END=1, err=1) uwork(1:10)
253 n_urpoly = 0
254 IF (dot_product(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
255 n_urpoly = 1
256 DO k = 2, 9
257 IF (abs(uwork(k)) >= 1.e-12_dp) n_urpoly = k
258 END DO
259 END IF
260! Polynomials of length 1 are not allowed, it seems we should use spline after all
261! This is creative guessing!
262 IF (n_urpoly < 2) n_urpoly = 0
263 END IF
264
265 CALL para_env%bcast(n_urpoly)
266 CALL para_env%bcast(uwork)
267 CALL para_env%bcast(ngrd)
268 CALL para_env%bcast(dgrd)
269
270 CALL para_env%bcast(skself)
271 CALL para_env%bcast(energy)
272 CALL para_env%bcast(eta)
273 CALL para_env%bcast(occupation)
274
275 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, &
276 z=z, zeff=sum(occupation), defined=.true., &
277 skself=skself, energy=energy, eta=eta, occupation=occupation)
278
279 ! Slater-Koster table
280 ALLOCATE (fmat(ngrd, 10))
281 ALLOCATE (smat(ngrd, 10))
282 IF (para_env%is_source()) THEN
283 DO k = 1, ngrd
284 READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
285 fmat(k, 1:10) = fwork(1:10)
286 smat(k, 1:10) = swork(1:10)
287 END DO
288 END IF
289 CALL para_env%bcast(fmat)
290 CALL para_env%bcast(smat)
291
292 !
293 ! Determine lmax for atom type.
294 ! An atomic orbital is 'active' if either its onsite energy is different from zero,
295 ! or
296 ! if this matrix element contains non-zero elements.
297 ! The sigma interactions are sufficient for that.
298 ! In the DFTB-Slako convention they are on orbital 10 (s-s-sigma),
299 ! 7 (p-p-sigma) and 3 (d-d-sigma).
300 !
301 ! We also allow lmax to be set in the input (in KIND)
302 !
303 CALL get_qs_kind(qs_kind_set(ikind), lmax_dftb=lmax)
304 IF (lmax < 0) THEN
305 lmax = 0
306 DO l = 0, 3
307 SELECT CASE (l)
308 CASE DEFAULT
309 cpabort("")
310 CASE (0)
311 lp = 10
312 CASE (1)
313 lp = 7
314 CASE (2)
315 lp = 3
316 CASE (3)
317 lp = 3 ! this is wrong but we don't allow f anyway
318 END SELECT
319 ! Technical note: In some slako files dummies are included in the
320 ! first matrix elements, so remove them.
321 IF ((abs(skself(l)) > 0._dp) .OR. &
322 (sum(abs(fmat(ngrd/10:ngrd, lp))) > 0._dp)) lmax = l
323 END DO
324 ! l=2 (d) is maximum
325 lmax = min(2, lmax)
326 END IF
327 IF (lmax > 2) THEN
328 CALL cp_abort(__location__, "Maximum L allowed is d. "// &
329 "Use KIND/LMAX_DFTB to set smaller values if needed.")
330 END IF
331 !
332 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, &
333 lmax=lmax, natorb=(lmax + 1)**2)
334
335 spdim = 0
336 IF (n_urpoly == 0) THEN
337 IF (para_env%is_source()) THEN
338 ! Look for spline representation of repulsive potential
339 search = .true.
340 DO WHILE (search)
341 READ (runit, fmt='(A6)', END=1, err=1) cspline
342 IF (cspline == 'Spline') THEN
343 search = .false.
344 ! spline dimension and left-hand cutoff
345 READ (runit, fmt=*, END=1, err=1) spdim, s_cut
346 ALLOCATE (spxr(spdim, 2))
347 ALLOCATE (scoeff(spdim, 4))
348 ! e-functions describing left-hand extrapolation
349 READ (runit, fmt=*, END=1, err=1) srep(1:3)
350 DO isp = 1, spdim - 1
351 ! location and coefficients of 'normal' spline range
352 READ (runit, fmt=*, END=1, err=1) spxr(isp, 1:2), scoeff(isp, 1:4)
353 END DO
354 ! last point has 2 more coefficients
355 READ (runit, fmt=*, END=1, err=1) spxr(spdim, 1:2), scoeff(spdim, 1:4), surr(1:2)
356 END IF
357 END DO
358 END IF
359 END IF
360
361 IF (para_env%is_source()) THEN
362 CALL close_file(unit_number=runit)
363 END IF
364
365 CALL para_env%bcast(spdim)
366 IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
367 ALLOCATE (spxr(spdim, 2))
368 ALLOCATE (scoeff(spdim, 4))
369 END IF
370 IF (spdim > 0) THEN
371 CALL para_env%bcast(spxr)
372 CALL para_env%bcast(scoeff)
373 CALL para_env%bcast(surr)
374 CALL para_env%bcast(srep)
375 CALL para_env%bcast(s_cut)
376 END IF
377
378 ! store potential data
379 ! allocate data
380 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
381 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
382 llm = 0
383 DO l1 = 0, max(lmax_a, lmax_b)
384 DO l2 = 0, min(l1, lmax_a, lmax_b)
385 DO m = 0, l2
386 llm = llm + 1
387 END DO
388 END DO
389 END DO
390 CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
391 ngrd, llm, spdim)
392
393 ! repulsive potential
394 dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
395 dftb_potential(ikind, jkind)%urep_cut = uwork(10)
396 dftb_potential(ikind, jkind)%urep(:) = 0._dp
397 dftb_potential(ikind, jkind)%urep(1) = uwork(10)
398 dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
399
400 ! Slater-Koster tables
401 dftb_potential(ikind, jkind)%dgrd = dgrd
402 CALL skreorder(fmat, lmax_a, lmax_b)
403 dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
404 CALL skreorder(smat, lmax_a, lmax_b)
405 dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
406 dftb_potential(ikind, jkind)%ngrdcut = ngrd + int(slako_d0/dgrd)
407 ! Splines
408 IF (spdim > 0) THEN
409 dftb_potential(ikind, jkind)%s_cut = s_cut
410 dftb_potential(ikind, jkind)%srep = srep
411 dftb_potential(ikind, jkind)%spxr = spxr
412 dftb_potential(ikind, jkind)%scoeff = scoeff
413 dftb_potential(ikind, jkind)%surr = surr
414 END IF
415
416 DEALLOCATE (fmat)
417 DEALLOCATE (smat)
418 IF (spdim > 0) THEN
419 DEALLOCATE (spxr)
420 DEALLOCATE (scoeff)
421 END IF
422
423 END DO
424
425 ! no all other pairs
426 DO ikind = 1, nkind
427 CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
428 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
429
430 IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
431 CALL allocate_dftb_atom_param(dftb_atom_a)
432 CALL set_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
433 END IF
434
435 DO jkind = 1, nkind
436
437 IF (ikind == jkind) cycle
438 CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
439 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
440
441 IF (output_unit > 0) THEN
442 WRITE (output_unit, "(A,T30,A50)") " DFTB| Reading parameter file ", &
443 adjustr(trim(sk_files(ikind, jkind)))
444 END IF
445 skself = 0._dp
446 eta = 0._dp
447 occupation = 0._dp
448 IF (para_env%is_source()) THEN
449 runit = get_unit_number()
450 CALL open_file(file_name=sk_files(ikind, jkind), unit_number=runit)
451 ! grid density and number of grid poin ts
452 READ (runit, fmt=*, END=1, err=1) dgrd, ngrd
453!
454! ngrd -1 ?
455! In Slako tables, the grid starts at 0.0, in deMon it starts with dgrd
456!
457 ngrd = ngrd - 1
458!
459 IF (ikind == jkind) THEN
460 ! orbital energy, total energy, hardness, occupation
461 READ (runit, fmt=*, END=1, err=1) skself(2:0:-1), energy, &
462 eta(2:0:-1), occupation(2:0:-1)
463 END IF
464 ! repulsive potential as polynomial
465 READ (runit, fmt=*, END=1, err=1) uwork(1:10)
466 n_urpoly = 0
467 IF (dot_product(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
468 n_urpoly = 1
469 DO k = 2, 9
470 IF (abs(uwork(k)) >= 1.e-12_dp) n_urpoly = k
471 END DO
472 END IF
473! Polynomials of length 1 are not allowed, it seems we should use spline after all
474! This is creative guessing!
475 IF (n_urpoly < 2) n_urpoly = 0
476 END IF
477
478 CALL para_env%bcast(n_urpoly)
479 CALL para_env%bcast(uwork)
480 CALL para_env%bcast(ngrd)
481 CALL para_env%bcast(dgrd)
482
483 ! Slater-Koster table
484 ALLOCATE (fmat(ngrd, 10))
485 ALLOCATE (smat(ngrd, 10))
486 IF (para_env%is_source()) THEN
487 DO k = 1, ngrd
488 READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
489 fmat(k, 1:10) = fwork(1:10)
490 smat(k, 1:10) = swork(1:10)
491 END DO
492 END IF
493 CALL para_env%bcast(fmat)
494 CALL para_env%bcast(smat)
495
496 spdim = 0
497 IF (n_urpoly == 0) THEN
498 IF (para_env%is_source()) THEN
499 ! Look for spline representation of repulsive potential
500 search = .true.
501 DO WHILE (search)
502 READ (runit, fmt='(A6)', END=1, err=1) cspline
503 IF (cspline == 'Spline') THEN
504 search = .false.
505 ! spline dimension and left-hand cutoff
506 READ (runit, fmt=*, END=1, err=1) spdim, s_cut
507 ALLOCATE (spxr(spdim, 2))
508 ALLOCATE (scoeff(spdim, 4))
509 ! e-functions describing left-hand extrapolation
510 READ (runit, fmt=*, END=1, err=1) srep(1:3)
511 DO isp = 1, spdim - 1
512 ! location and coefficients of 'normal' spline range
513 READ (runit, fmt=*, END=1, err=1) spxr(isp, 1:2), scoeff(isp, 1:4)
514 END DO
515 ! last point has 2 more coefficients
516 READ (runit, fmt=*, END=1, err=1) spxr(spdim, 1:2), scoeff(spdim, 1:4), surr(1:2)
517 END IF
518 END DO
519 END IF
520 END IF
521
522 IF (para_env%is_source()) THEN
523 CALL close_file(unit_number=runit)
524 END IF
525
526 CALL para_env%bcast(spdim)
527 IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
528 ALLOCATE (spxr(spdim, 2))
529 ALLOCATE (scoeff(spdim, 4))
530 END IF
531 IF (spdim > 0) THEN
532 CALL para_env%bcast(spxr)
533 CALL para_env%bcast(scoeff)
534 CALL para_env%bcast(surr)
535 CALL para_env%bcast(srep)
536 CALL para_env%bcast(s_cut)
537 END IF
538
539 ! store potential data
540 ! allocate data
541 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
542 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
543 llm = 0
544 DO l1 = 0, max(lmax_a, lmax_b)
545 DO l2 = 0, min(l1, lmax_a, lmax_b)
546 DO m = 0, l2
547 llm = llm + 1
548 END DO
549 END DO
550 END DO
551 CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
552 ngrd, llm, spdim)
553
554 ! repulsive potential
555 dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
556 dftb_potential(ikind, jkind)%urep_cut = uwork(10)
557 dftb_potential(ikind, jkind)%urep(:) = 0._dp
558 dftb_potential(ikind, jkind)%urep(1) = uwork(10)
559 dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
560
561 ! Slater-Koster tables
562 dftb_potential(ikind, jkind)%dgrd = dgrd
563 CALL skreorder(fmat, lmax_a, lmax_b)
564 dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
565 CALL skreorder(smat, lmax_a, lmax_b)
566 dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
567 dftb_potential(ikind, jkind)%ngrdcut = ngrd + int(slako_d0/dgrd)
568 ! Splines
569 IF (spdim > 0) THEN
570 dftb_potential(ikind, jkind)%s_cut = s_cut
571 dftb_potential(ikind, jkind)%srep = srep
572 dftb_potential(ikind, jkind)%spxr = spxr
573 dftb_potential(ikind, jkind)%scoeff = scoeff
574 dftb_potential(ikind, jkind)%surr = surr
575 END IF
576
577 DEALLOCATE (fmat)
578 DEALLOCATE (smat)
579 IF (spdim > 0) THEN
580 DEALLOCATE (spxr)
581 DEALLOCATE (scoeff)
582 END IF
583
584 END DO
585 END DO
586
587 DEALLOCATE (sk_files)
588
589 ! read dispersion parameters (UFF type)
590 IF (dftb_control%dispersion) THEN
591
592 IF (dftb_control%dispersion_type == dispersion_uff) THEN
593 file_name = trim(dftb_control%sk_file_path)//"/"// &
594 trim(dftb_control%uff_force_field)
595 block
596 TYPE(cp_parser_type) :: parser
597 DO ikind = 1, nkind
598 CALL get_atomic_kind(atomic_kind_set(ikind), name=iname)
599 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
600
601 m = len_trim(iname)
602 CALL parser_create(parser, file_name, para_env=para_env)
603 found = .false.
604 DO
605 at_end = .false.
606 CALL parser_get_next_line(parser, 1, at_end)
607 IF (at_end) EXIT
608 CALL parser_get_object(parser, name_a)
609 ! parser is no longer removing leading quotes
610 IF (name_a(1:1) == '"') name_a(1:m) = name_a(2:m + 1)
611 IF (name_a(1:m) == trim(iname)) THEN
612 CALL parser_get_object(parser, rb)
613 CALL parser_get_object(parser, rb)
614 CALL parser_get_object(parser, ra)
615 CALL parser_get_object(parser, da)
616 found = .true.
617 ra = ra/angstrom
618 da = da/kcalmol
619 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, name=iname, xi=ra, di=da)
620 EXIT
621 END IF
622 END DO
623 CALL parser_release(parser)
624 END DO
625 END block
626 END IF
627
628 END IF
629
630 ! extract simple atom interaction radii
631 DO ikind = 1, nkind
632 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
633 radmax = (dftb_potential(ikind, ikind)%ngrdcut + 1)* &
634 dftb_potential(ikind, ikind)%dgrd*0.5_dp
635 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=radmax)
636 END DO
637 DO ikind = 1, nkind
638 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
639 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
640 DO jkind = 1, nkind
641 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
642 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, cutoff=rb)
643 radmax = (dftb_potential(ikind, jkind)%ngrdcut + 1)* &
644 dftb_potential(ikind, jkind)%dgrd
645 IF (ra + rb < radmax) THEN
646 ra = ra + (radmax - ra - rb)*0.5_dp
647 rb = rb + (radmax - ra - rb)*0.5_dp
648 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
649 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_b, cutoff=rb)
650 END IF
651 END DO
652 END DO
653
654 ! set correct core charge in potential
655 DO ikind = 1, nkind
656 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
657 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, zeff=zeff)
658 CALL set_potential(potential=qs_kind_set(ikind)%all_potential, &
659 zeff=zeff, zeff_correction=0.0_dp)
660 END DO
661
662 ! setup DFTB3 parameters
663 IF (dftb_control%dftb3_diagonal) THEN
664 DO ikind = 1, nkind
665 CALL get_qs_kind(qs_kind_set(ikind), dftb3_param=db)
666 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
667 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, dudq=db)
668 END DO
669 END IF
670
671 ! setup dispersion parameters (UFF type)
672 IF (dftb_control%dispersion) THEN
673 IF (dftb_control%dispersion_type == dispersion_uff) THEN
674 eps_disp = dftb_control%eps_disp
675 DO ikind = 1, nkind
676 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
677 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, xi=ra, di=da)
678 rcdisp = 0._dp
679 DO jkind = 1, nkind
680 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
681 CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, xi=rb, di=db)
682 xij = sqrt(ra*rb)
683 dij = sqrt(da*db)
684 dftb_potential(ikind, jkind)%xij = xij
685 dftb_potential(ikind, jkind)%dij = dij
686 dftb_potential(ikind, jkind)%x0ij = xij*(0.5_dp**(1.0_dp/6.0_dp))
687 dftb_potential(ikind, jkind)%a = dij*396.0_dp/25.0_dp
688 dftb_potential(ikind, jkind)%b = &
689 dij/(xij**5)*672.0_dp*2.0_dp**(5.0_dp/6.0_dp)/25.0_dp
690 dftb_potential(ikind, jkind)%c = &
691 -dij/(xij**10)*2.0_dp**(2.0_dp/3.0_dp)*552.0_dp/25.0_dp
692 rmax6 = ((8._dp*pi*dij/eps_disp)*xij**6)**0.25_dp
693 rcdisp = max(rcdisp, rmax6*0.5_dp)
694 END DO
695 CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, rcdisp=rcdisp)
696 END DO
697 END IF
698 END IF
699
700 RETURN
701
7021 CONTINUE
703 cpabort("")
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:296
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:308
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:237
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F: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:1063
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:17
input_constants::dispersion_uff
integer, parameter, public dispersion_uff
Definition input_constants.F:1207
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:177
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:123
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, 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, 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:454
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:2840
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:111
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:763
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:171