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