Line data Source code
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 : ! **************************************************************************************************
11 : MODULE 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 :
68 : CONTAINS
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 292 : 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 292 : 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 292 : 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 292 : output_unit = -1
112 292 : NULLIFY (logger)
113 292 : logger => cp_get_default_logger()
114 292 : 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 0 : "PRINT%KINDS", extension=".Log")
118 0 : IF (output_unit > 0) THEN
119 : WRITE (output_unit, "(/,A)") " DFTB| A set of relativistic DFTB "// &
120 0 : "parameters for material sciences."
121 : WRITE (output_unit, "(A)") " DFTB| J. Frenzel, N. Jardillier, A.F. Oliveira,"// &
122 0 : " T. Heine, G. Seifert"
123 0 : WRITE (output_unit, "(A)") " DFTB| TU Dresden, 2002-2007"
124 0 : WRITE (output_unit, "(/,A)") " DFTB| Non-SCC parameters "
125 0 : WRITE (output_unit, "(A,T25,A)") " DFTB| C,H :", &
126 0 : " D. Porezag et al, PRB 51 12947 (1995)"
127 0 : WRITE (output_unit, "(A,T25,A)") " DFTB| B,N :", &
128 0 : " J. Widany et al, PRB 53 4443 (1996)"
129 0 : WRITE (output_unit, "(A,T25,A)") " DFTB| Li,Na,K,Cl :", &
130 0 : " S. Hazebroucq et al, JCP 123 134510 (2005)"
131 0 : WRITE (output_unit, "(A,T25,A)") " DFTB| F :", &
132 0 : " T. Heine et al, JCSoc-Perkins Trans 2 707 (1999)"
133 0 : WRITE (output_unit, "(A,T25,A)") " DFTB| Mo,S :", &
134 0 : " G. Seifert et al, PRL 85 146 (2000)"
135 0 : WRITE (output_unit, "(A,T25,A)") " DFTB| P :", &
136 0 : " G. Seifert et al, EPS 16 341 (2001)"
137 0 : WRITE (output_unit, "(A,T25,A)") " DFTB| Sc,N,C :", &
138 0 : " M. Krause et al, JCP 115 6596 (2001)"
139 : END IF
140 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
141 0 : "PRINT%KINDS")
142 : END IF
143 :
144 292 : sklist = (dftb_control%sk_file_list /= "")
145 :
146 292 : nkind = SIZE(atomic_kind_set)
147 1168 : ALLOCATE (sk_files(nkind, nkind))
148 : ! allocate potential structures
149 7588 : ALLOCATE (dftb_potential(nkind, nkind))
150 292 : CALL qs_dftb_pairpot_init(dftb_potential)
151 :
152 910 : DO ikind = 1, nkind
153 618 : CALL get_atomic_kind(atomic_kind_set(ikind), name=iname, element_symbol=iel)
154 618 : CALL uppercase(iname)
155 618 : CALL uppercase(iel)
156 618 : ldum = qmmm_ff_precond_only_qm(iname)
157 2332 : DO jkind = 1, nkind
158 1422 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname, element_symbol=jel)
159 1422 : CALL uppercase(jname)
160 1422 : CALL uppercase(jel)
161 1422 : ldum = qmmm_ff_precond_only_qm(jname)
162 1422 : found = .FALSE.
163 1518 : DO k = 1, SIZE(dftb_control%sk_pair_list, 2)
164 192 : name_a = TRIM(dftb_control%sk_pair_list(1, k))
165 192 : name_b = TRIM(dftb_control%sk_pair_list(2, k))
166 192 : CALL uppercase(name_a)
167 192 : CALL uppercase(name_b)
168 1518 : IF ((iname == name_a .AND. jname == name_b)) THEN
169 : sk_files(ikind, jkind) = TRIM(dftb_control%sk_file_path)//"/"// &
170 96 : TRIM(dftb_control%sk_pair_list(3, k))
171 96 : found = .TRUE.
172 96 : EXIT
173 : END IF
174 : END DO
175 1422 : IF (.NOT. found .AND. sklist) THEN
176 : file_name = TRIM(dftb_control%sk_file_path)//"/"// &
177 1326 : TRIM(dftb_control%sk_file_list)
178 1326 : BLOCK
179 : TYPE(cp_parser_type) :: parser
180 1326 : CALL parser_create(parser, file_name, para_env=para_env)
181 : DO
182 : at_end = .FALSE.
183 20066 : CALL parser_get_next_line(parser, 1, at_end)
184 20066 : IF (at_end) EXIT
185 20066 : CALL parser_get_object(parser, name_a, lower_to_upper=.TRUE.)
186 20066 : CALL parser_get_object(parser, name_b, lower_to_upper=.TRUE.)
187 : !Checking Names
188 20066 : IF ((iname == name_a .AND. jname == name_b)) THEN
189 1326 : CALL parser_get_object(parser, skfn, string_length=8)
190 : sk_files(ikind, jkind) = TRIM(dftb_control%sk_file_path)//"/"// &
191 1326 : TRIM(skfn)
192 1326 : found = .TRUE.
193 1326 : EXIT
194 : END IF
195 : !Checking Element
196 18740 : IF ((iel == name_a .AND. jel == name_b)) THEN
197 0 : CALL parser_get_object(parser, skfn, string_length=8)
198 : sk_files(ikind, jkind) = TRIM(dftb_control%sk_file_path)//"/"// &
199 0 : TRIM(skfn)
200 0 : found = .TRUE.
201 0 : EXIT
202 : END IF
203 : END DO
204 5304 : CALL parser_release(parser)
205 : END BLOCK
206 : END IF
207 2040 : IF (.NOT. found) THEN
208 : CALL cp_abort(__LOCATION__, &
209 : "Failure in assigning KINDS <"//TRIM(iname)//"> and <"//TRIM(jname)// &
210 0 : "> 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 910 : DO ikind = 1, nkind
217 618 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
218 :
219 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
220 618 : IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
221 618 : CALL allocate_dftb_atom_param(dftb_atom_a)
222 618 : 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 618 : jkind = ikind
227 :
228 618 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
229 618 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
230 :
231 618 : IF (output_unit > 0) THEN
232 0 : WRITE (output_unit, "(A,T30,A50)") " DFTB| Reading parameter file ", &
233 0 : ADJUSTR(TRIM(sk_files(jkind, ikind)))
234 : END IF
235 618 : skself = 0._dp
236 618 : eta = 0._dp
237 618 : occupation = 0._dp
238 618 : IF (para_env%is_source()) THEN
239 309 : runit = get_unit_number()
240 309 : CALL open_file(file_name=sk_files(jkind, ikind), unit_number=runit)
241 : ! grid density and number of grid poin ts
242 309 : 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 309 : ngrd = ngrd - 1
248 : !
249 : ! orbital energy, total energy, hardness, occupation
250 309 : READ (runit, fmt=*, END=1, err=1) skself(2:0:-1), energy, &
251 618 : eta(2:0:-1), occupation(2:0:-1)
252 : ! repulsive potential as polynomial
253 309 : READ (runit, fmt=*, END=1, err=1) uwork(1:10)
254 309 : n_urpoly = 0
255 3090 : IF (DOT_PRODUCT(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
256 71 : n_urpoly = 1
257 639 : DO k = 2, 9
258 639 : 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 309 : IF (n_urpoly < 2) n_urpoly = 0
264 : END IF
265 :
266 618 : CALL para_env%bcast(n_urpoly)
267 618 : CALL para_env%bcast(uwork)
268 618 : CALL para_env%bcast(ngrd)
269 618 : CALL para_env%bcast(dgrd)
270 :
271 618 : CALL para_env%bcast(skself)
272 618 : CALL para_env%bcast(energy)
273 618 : CALL para_env%bcast(eta)
274 618 : 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 3090 : skself=skself, energy=energy, eta=eta, occupation=occupation)
279 :
280 : ! Slater-Koster table
281 1854 : ALLOCATE (fmat(ngrd, 10))
282 1236 : ALLOCATE (smat(ngrd, 10))
283 618 : IF (para_env%is_source()) THEN
284 152900 : DO k = 1, ngrd
285 152591 : READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
286 1678501 : fmat(k, 1:10) = fwork(1:10)
287 1678810 : smat(k, 1:10) = swork(1:10)
288 : END DO
289 : END IF
290 618 : CALL para_env%bcast(fmat)
291 618 : 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 618 : CALL get_qs_kind(qs_kind_set(ikind), lmax_dftb=lmax)
305 618 : IF (lmax < 0) THEN
306 614 : lmax = 0
307 3070 : DO l = 0, 3
308 0 : SELECT CASE (l)
309 : CASE DEFAULT
310 0 : CPABORT("Only 0, 1, 2 are supported as the value of l")
311 : CASE (0)
312 614 : lp = 10
313 : CASE (1)
314 614 : lp = 7
315 : CASE (2)
316 614 : lp = 3
317 : CASE (3)
318 2456 : 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 1098592 : IF ((ABS(skself(l)) > 0._dp) .OR. &
323 1646 : (SUM(ABS(fmat(ngrd/10:ngrd, lp))) > 0._dp)) lmax = l
324 : END DO
325 : ! l=2 (d) is maximum
326 614 : lmax = MIN(2, lmax)
327 : END IF
328 618 : IF (lmax > 2) THEN
329 : CALL cp_abort(__LOCATION__, "Maximum L allowed is d. "// &
330 0 : "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 618 : lmax=lmax, natorb=(lmax + 1)**2)
335 :
336 618 : spdim = 0
337 618 : IF (n_urpoly == 0) THEN
338 476 : IF (para_env%is_source()) THEN
339 : ! Look for spline representation of repulsive potential
340 : search = .TRUE.
341 : DO WHILE (search)
342 5000 : READ (runit, fmt='(A6)', END=1, err=1) cspline
343 5000 : IF (cspline == 'Spline') THEN
344 238 : search = .FALSE.
345 : ! spline dimension and left-hand cutoff
346 238 : READ (runit, fmt=*, END=1, err=1) spdim, s_cut
347 714 : ALLOCATE (spxr(spdim, 2))
348 714 : ALLOCATE (scoeff(spdim, 4))
349 : ! e-functions describing left-hand extrapolation
350 238 : READ (runit, fmt=*, END=1, err=1) srep(1:3)
351 7722 : DO isp = 1, spdim - 1
352 : ! location and coefficients of 'normal' spline range
353 7722 : 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 238 : 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 618 : IF (para_env%is_source()) THEN
363 309 : CALL close_file(unit_number=runit)
364 : END IF
365 :
366 618 : CALL para_env%bcast(spdim)
367 618 : IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
368 714 : ALLOCATE (spxr(spdim, 2))
369 714 : ALLOCATE (scoeff(spdim, 4))
370 : END IF
371 618 : IF (spdim > 0) THEN
372 476 : CALL para_env%bcast(spxr)
373 476 : CALL para_env%bcast(scoeff)
374 476 : CALL para_env%bcast(surr)
375 476 : CALL para_env%bcast(srep)
376 476 : CALL para_env%bcast(s_cut)
377 : END IF
378 :
379 : ! store potential data
380 : ! allocate data
381 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
382 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
383 618 : llm = 0
384 1634 : DO l1 = 0, MAX(lmax_a, lmax_b)
385 3072 : DO l2 = 0, MIN(l1, lmax_a, lmax_b)
386 4338 : DO m = 0, l2
387 3322 : llm = llm + 1
388 : END DO
389 : END DO
390 : END DO
391 : CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
392 618 : ngrd, llm, spdim)
393 :
394 : ! repulsive potential
395 618 : dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
396 618 : dftb_potential(ikind, jkind)%urep_cut = uwork(10)
397 6798 : dftb_potential(ikind, jkind)%urep(:) = 0._dp
398 618 : dftb_potential(ikind, jkind)%urep(1) = uwork(10)
399 1526 : dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
400 :
401 : ! Slater-Koster tables
402 618 : dftb_potential(ikind, jkind)%dgrd = dgrd
403 618 : CALL skreorder(fmat, lmax_a, lmax_b)
404 923578 : dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
405 618 : CALL skreorder(smat, lmax_a, lmax_b)
406 923578 : dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
407 618 : dftb_potential(ikind, jkind)%ngrdcut = ngrd + INT(slako_d0/dgrd)
408 : ! Splines
409 618 : IF (spdim > 0) THEN
410 476 : dftb_potential(ikind, jkind)%s_cut = s_cut
411 1904 : dftb_potential(ikind, jkind)%srep = srep
412 32316 : dftb_potential(ikind, jkind)%spxr = spxr
413 64156 : dftb_potential(ikind, jkind)%scoeff = scoeff
414 1428 : dftb_potential(ikind, jkind)%surr = surr
415 : END IF
416 :
417 618 : DEALLOCATE (fmat)
418 618 : DEALLOCATE (smat)
419 2764 : IF (spdim > 0) THEN
420 476 : DEALLOCATE (spxr)
421 476 : DEALLOCATE (scoeff)
422 : END IF
423 :
424 : END DO
425 :
426 : ! no all other pairs
427 910 : DO ikind = 1, nkind
428 618 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
429 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
430 :
431 618 : IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
432 0 : CALL allocate_dftb_atom_param(dftb_atom_a)
433 0 : CALL set_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
434 : END IF
435 :
436 2950 : DO jkind = 1, nkind
437 :
438 1422 : IF (ikind == jkind) CYCLE
439 804 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
440 804 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
441 :
442 804 : IF (output_unit > 0) THEN
443 0 : WRITE (output_unit, "(A,T30,A50)") " DFTB| Reading parameter file ", &
444 0 : ADJUSTR(TRIM(sk_files(ikind, jkind)))
445 : END IF
446 804 : skself = 0._dp
447 804 : eta = 0._dp
448 804 : occupation = 0._dp
449 804 : IF (para_env%is_source()) THEN
450 402 : runit = get_unit_number()
451 402 : CALL open_file(file_name=sk_files(ikind, jkind), unit_number=runit)
452 : ! grid density and number of grid poin ts
453 402 : 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 402 : 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 402 : READ (runit, fmt=*, END=1, err=1) uwork(1:10)
467 402 : n_urpoly = 0
468 4020 : IF (DOT_PRODUCT(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
469 86 : n_urpoly = 1
470 774 : DO k = 2, 9
471 774 : 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 402 : IF (n_urpoly < 2) n_urpoly = 0
477 : END IF
478 :
479 804 : CALL para_env%bcast(n_urpoly)
480 804 : CALL para_env%bcast(uwork)
481 804 : CALL para_env%bcast(ngrd)
482 804 : CALL para_env%bcast(dgrd)
483 :
484 : ! Slater-Koster table
485 2412 : ALLOCATE (fmat(ngrd, 10))
486 1608 : ALLOCATE (smat(ngrd, 10))
487 804 : IF (para_env%is_source()) THEN
488 199960 : DO k = 1, ngrd
489 199558 : READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
490 2195138 : fmat(k, 1:10) = fwork(1:10)
491 2195540 : smat(k, 1:10) = swork(1:10)
492 : END DO
493 : END IF
494 804 : CALL para_env%bcast(fmat)
495 804 : CALL para_env%bcast(smat)
496 :
497 804 : spdim = 0
498 804 : IF (n_urpoly == 0) THEN
499 632 : IF (para_env%is_source()) THEN
500 : ! Look for spline representation of repulsive potential
501 : search = .TRUE.
502 : DO WHILE (search)
503 6444 : READ (runit, fmt='(A6)', END=1, err=1) cspline
504 6444 : IF (cspline == 'Spline') THEN
505 316 : search = .FALSE.
506 : ! spline dimension and left-hand cutoff
507 316 : READ (runit, fmt=*, END=1, err=1) spdim, s_cut
508 948 : ALLOCATE (spxr(spdim, 2))
509 948 : ALLOCATE (scoeff(spdim, 4))
510 : ! e-functions describing left-hand extrapolation
511 316 : READ (runit, fmt=*, END=1, err=1) srep(1:3)
512 10156 : DO isp = 1, spdim - 1
513 : ! location and coefficients of 'normal' spline range
514 10156 : 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 316 : 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 804 : IF (para_env%is_source()) THEN
524 402 : CALL close_file(unit_number=runit)
525 : END IF
526 :
527 804 : CALL para_env%bcast(spdim)
528 804 : IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
529 948 : ALLOCATE (spxr(spdim, 2))
530 948 : ALLOCATE (scoeff(spdim, 4))
531 : END IF
532 804 : IF (spdim > 0) THEN
533 632 : CALL para_env%bcast(spxr)
534 632 : CALL para_env%bcast(scoeff)
535 632 : CALL para_env%bcast(surr)
536 632 : CALL para_env%bcast(srep)
537 632 : CALL para_env%bcast(s_cut)
538 : END IF
539 :
540 : ! store potential data
541 : ! allocate data
542 804 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
543 804 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
544 804 : llm = 0
545 2428 : DO l1 = 0, MAX(lmax_a, lmax_b)
546 4260 : DO l2 = 0, MIN(l1, lmax_a, lmax_b)
547 5504 : DO m = 0, l2
548 3880 : llm = llm + 1
549 : END DO
550 : END DO
551 : END DO
552 : CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
553 804 : ngrd, llm, spdim)
554 :
555 : ! repulsive potential
556 804 : dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
557 804 : dftb_potential(ikind, jkind)%urep_cut = uwork(10)
558 8844 : dftb_potential(ikind, jkind)%urep(:) = 0._dp
559 804 : dftb_potential(ikind, jkind)%urep(1) = uwork(10)
560 1748 : dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
561 :
562 : ! Slater-Koster tables
563 804 : dftb_potential(ikind, jkind)%dgrd = dgrd
564 804 : CALL skreorder(fmat, lmax_a, lmax_b)
565 1012644 : dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
566 804 : CALL skreorder(smat, lmax_a, lmax_b)
567 1012644 : dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
568 804 : dftb_potential(ikind, jkind)%ngrdcut = ngrd + INT(slako_d0/dgrd)
569 : ! Splines
570 804 : IF (spdim > 0) THEN
571 632 : dftb_potential(ikind, jkind)%s_cut = s_cut
572 2528 : dftb_potential(ikind, jkind)%srep = srep
573 42520 : dftb_potential(ikind, jkind)%spxr = spxr
574 84408 : dftb_potential(ikind, jkind)%scoeff = scoeff
575 1896 : dftb_potential(ikind, jkind)%surr = surr
576 : END IF
577 :
578 804 : DEALLOCATE (fmat)
579 804 : DEALLOCATE (smat)
580 2226 : IF (spdim > 0) THEN
581 632 : DEALLOCATE (spxr)
582 632 : DEALLOCATE (scoeff)
583 : END IF
584 :
585 : END DO
586 : END DO
587 :
588 292 : DEALLOCATE (sk_files)
589 :
590 : ! read dispersion parameters (UFF type)
591 292 : IF (dftb_control%dispersion) THEN
592 :
593 112 : IF (dftb_control%dispersion_type == dispersion_uff) THEN
594 : file_name = TRIM(dftb_control%sk_file_path)//"/"// &
595 94 : TRIM(dftb_control%uff_force_field)
596 : BLOCK
597 : TYPE(cp_parser_type) :: parser
598 476 : DO ikind = 1, nkind
599 194 : CALL get_atomic_kind(atomic_kind_set(ikind), name=iname)
600 194 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
601 :
602 194 : m = LEN_TRIM(iname)
603 194 : CALL parser_create(parser, file_name, para_env=para_env)
604 194 : found = .FALSE.
605 : DO
606 : at_end = .FALSE.
607 2104 : CALL parser_get_next_line(parser, 1, at_end)
608 2104 : IF (at_end) EXIT
609 2104 : CALL parser_get_object(parser, name_a)
610 : ! parser is no longer removing leading quotes
611 2104 : IF (name_a(1:1) == '"') name_a(1:m) = name_a(2:m + 1)
612 2104 : IF (name_a(1:m) == TRIM(iname)) THEN
613 194 : CALL parser_get_object(parser, rb)
614 194 : CALL parser_get_object(parser, rb)
615 194 : CALL parser_get_object(parser, ra)
616 194 : CALL parser_get_object(parser, da)
617 194 : found = .TRUE.
618 194 : ra = ra/angstrom
619 194 : da = da/kcalmol
620 194 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, name=iname, xi=ra, di=da)
621 194 : EXIT
622 : END IF
623 : END DO
624 482 : 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 910 : DO ikind = 1, nkind
633 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
634 : radmax = (dftb_potential(ikind, ikind)%ngrdcut + 1)* &
635 618 : dftb_potential(ikind, ikind)%dgrd*0.5_dp
636 910 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=radmax)
637 : END DO
638 910 : DO ikind = 1, nkind
639 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
640 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
641 2332 : DO jkind = 1, nkind
642 1422 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
643 1422 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, cutoff=rb)
644 : radmax = (dftb_potential(ikind, jkind)%ngrdcut + 1)* &
645 1422 : dftb_potential(ikind, jkind)%dgrd
646 2040 : IF (ra + rb < radmax) THEN
647 8 : ra = ra + (radmax - ra - rb)*0.5_dp
648 8 : rb = rb + (radmax - ra - rb)*0.5_dp
649 8 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
650 8 : 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 910 : DO ikind = 1, nkind
657 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
658 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, zeff=zeff)
659 : CALL set_potential(potential=qs_kind_set(ikind)%all_potential, &
660 910 : zeff=zeff, zeff_correction=0.0_dp)
661 : END DO
662 :
663 : ! setup DFTB3 parameters
664 292 : IF (dftb_control%dftb3_diagonal) THEN
665 134 : DO ikind = 1, nkind
666 88 : CALL get_qs_kind(qs_kind_set(ikind), dftb3_param=db)
667 88 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
668 222 : 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 292 : IF (dftb_control%dispersion) THEN
674 112 : IF (dftb_control%dispersion_type == dispersion_uff) THEN
675 94 : eps_disp = dftb_control%eps_disp
676 288 : DO ikind = 1, nkind
677 194 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
678 194 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, xi=ra, di=da)
679 194 : rcdisp = 0._dp
680 652 : DO jkind = 1, nkind
681 458 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
682 458 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, xi=rb, di=db)
683 458 : xij = SQRT(ra*rb)
684 458 : dij = SQRT(da*db)
685 458 : dftb_potential(ikind, jkind)%xij = xij
686 458 : dftb_potential(ikind, jkind)%dij = dij
687 458 : dftb_potential(ikind, jkind)%x0ij = xij*(0.5_dp**(1.0_dp/6.0_dp))
688 458 : dftb_potential(ikind, jkind)%a = dij*396.0_dp/25.0_dp
689 : dftb_potential(ikind, jkind)%b = &
690 458 : 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 458 : -dij/(xij**10)*2.0_dp**(2.0_dp/3.0_dp)*552.0_dp/25.0_dp
693 458 : rmax6 = ((8._dp*pi*dij/eps_disp)*xij**6)**0.25_dp
694 652 : rcdisp = MAX(rcdisp, rmax6*0.5_dp)
695 : END DO
696 288 : 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 :
703 : 1 CONTINUE
704 : ! Many instances of READ (..., END=1, err=1) gets conflated here
705 : ! TODO: overhaul the file parser to handle errors separately
706 0 : CPABORT("Something went wrong while reading DFTB parameter file")
707 :
708 292 : 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 2844 : 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 1411440 : DO i = 1, SIZE(xmat, 1)
741 1408596 : skllm = 0._dp
742 1408596 : skllm(0, 0, 0) = xmat(i, 10)
743 1408596 : skllm(1, 0, 0) = xmat(i, 9)
744 1408596 : skllm(2, 0, 0) = xmat(i, 8)
745 1408596 : skllm(1, 1, 1) = xmat(i, 7)
746 1408596 : skllm(1, 1, 0) = xmat(i, 6)
747 1408596 : skllm(2, 1, 1) = xmat(i, 5)
748 1408596 : skllm(2, 1, 0) = xmat(i, 4)
749 1408596 : skllm(2, 2, 2) = xmat(i, 3)
750 1408596 : skllm(2, 2, 1) = xmat(i, 2)
751 1408596 : skllm(2, 2, 0) = xmat(i, 1)
752 1408596 : llm = 0
753 4018800 : DO l1 = 0, MAX(la, lb)
754 7234776 : DO l2 = 0, MIN(l1, la, lb)
755 9687916 : DO m = 0, l2
756 3861736 : llm = llm + 1
757 7080556 : xmat(i, llm) = skllm(l1, l2, m)
758 : END DO
759 : END DO
760 : END DO
761 : END DO
762 : !
763 2844 : END SUBROUTINE skreorder
764 :
765 : END MODULE qs_dftb_parameters
766 :
|
