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 1422 : IF (.NOT. found) &
208 : CALL cp_abort(__LOCATION__, &
209 : "Failure in assigning KINDS <"//TRIM(iname)//"> and <"//TRIM(jname)// &
210 618 : "> to a DFTB interaction pair!")
211 : END DO
212 : END DO
213 : ! reading the files
214 : ! read all pairs, equal kind first
215 910 : DO ikind = 1, nkind
216 618 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
217 :
218 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
219 618 : IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
220 618 : CALL allocate_dftb_atom_param(dftb_atom_a)
221 618 : 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 618 : jkind = ikind
226 :
227 618 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
228 618 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
229 :
230 618 : 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 618 : skself = 0._dp
235 618 : eta = 0._dp
236 618 : occupation = 0._dp
237 618 : IF (para_env%is_source()) THEN
238 309 : runit = get_unit_number()
239 309 : CALL open_file(file_name=sk_files(jkind, ikind), unit_number=runit)
240 : ! grid density and number of grid poin ts
241 309 : 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 309 : ngrd = ngrd - 1
247 : !
248 : ! orbital energy, total energy, hardness, occupation
249 309 : READ (runit, fmt=*, END=1, err=1) skself(2:0:-1), energy, &
250 618 : eta(2:0:-1), occupation(2:0:-1)
251 : ! repulsive potential as polynomial
252 309 : READ (runit, fmt=*, END=1, err=1) uwork(1:10)
253 309 : n_urpoly = 0
254 3090 : IF (DOT_PRODUCT(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
255 71 : n_urpoly = 1
256 639 : DO k = 2, 9
257 639 : 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 309 : IF (n_urpoly < 2) n_urpoly = 0
263 : END IF
264 :
265 618 : CALL para_env%bcast(n_urpoly)
266 618 : CALL para_env%bcast(uwork)
267 618 : CALL para_env%bcast(ngrd)
268 618 : CALL para_env%bcast(dgrd)
269 :
270 618 : CALL para_env%bcast(skself)
271 618 : CALL para_env%bcast(energy)
272 618 : CALL para_env%bcast(eta)
273 618 : 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 3090 : skself=skself, energy=energy, eta=eta, occupation=occupation)
278 :
279 : ! Slater-Koster table
280 1854 : ALLOCATE (fmat(ngrd, 10))
281 1236 : ALLOCATE (smat(ngrd, 10))
282 618 : IF (para_env%is_source()) THEN
283 152900 : DO k = 1, ngrd
284 152591 : READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
285 1678501 : fmat(k, 1:10) = fwork(1:10)
286 1678810 : smat(k, 1:10) = swork(1:10)
287 : END DO
288 : END IF
289 618 : CALL para_env%bcast(fmat)
290 618 : 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 618 : CALL get_qs_kind(qs_kind_set(ikind), lmax_dftb=lmax)
304 618 : IF (lmax < 0) THEN
305 614 : lmax = 0
306 3070 : DO l = 0, 3
307 0 : SELECT CASE (l)
308 : CASE DEFAULT
309 0 : CPABORT("")
310 : CASE (0)
311 614 : lp = 10
312 : CASE (1)
313 614 : lp = 7
314 : CASE (2)
315 614 : lp = 3
316 : CASE (3)
317 2456 : 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 1098592 : IF ((ABS(skself(l)) > 0._dp) .OR. &
322 1646 : (SUM(ABS(fmat(ngrd/10:ngrd, lp))) > 0._dp)) lmax = l
323 : END DO
324 : ! l=2 (d) is maximum
325 614 : lmax = MIN(2, lmax)
326 : END IF
327 618 : IF (lmax > 2) THEN
328 : CALL cp_abort(__LOCATION__, "Maximum L allowed is d. "// &
329 0 : "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 618 : lmax=lmax, natorb=(lmax + 1)**2)
334 :
335 618 : spdim = 0
336 618 : IF (n_urpoly == 0) THEN
337 476 : IF (para_env%is_source()) THEN
338 : ! Look for spline representation of repulsive potential
339 : search = .TRUE.
340 : DO WHILE (search)
341 5000 : READ (runit, fmt='(A6)', END=1, err=1) cspline
342 5000 : IF (cspline == 'Spline') THEN
343 238 : search = .FALSE.
344 : ! spline dimension and left-hand cutoff
345 238 : READ (runit, fmt=*, END=1, err=1) spdim, s_cut
346 714 : ALLOCATE (spxr(spdim, 2))
347 714 : ALLOCATE (scoeff(spdim, 4))
348 : ! e-functions describing left-hand extrapolation
349 238 : READ (runit, fmt=*, END=1, err=1) srep(1:3)
350 7722 : DO isp = 1, spdim - 1
351 : ! location and coefficients of 'normal' spline range
352 7722 : 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 238 : 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 618 : IF (para_env%is_source()) THEN
362 309 : CALL close_file(unit_number=runit)
363 : END IF
364 :
365 618 : CALL para_env%bcast(spdim)
366 618 : IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
367 714 : ALLOCATE (spxr(spdim, 2))
368 714 : ALLOCATE (scoeff(spdim, 4))
369 : END IF
370 618 : IF (spdim > 0) THEN
371 476 : CALL para_env%bcast(spxr)
372 476 : CALL para_env%bcast(scoeff)
373 476 : CALL para_env%bcast(surr)
374 476 : CALL para_env%bcast(srep)
375 476 : CALL para_env%bcast(s_cut)
376 : END IF
377 :
378 : ! store potential data
379 : ! allocate data
380 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
381 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
382 618 : llm = 0
383 1634 : DO l1 = 0, MAX(lmax_a, lmax_b)
384 3072 : DO l2 = 0, MIN(l1, lmax_a, lmax_b)
385 4338 : DO m = 0, l2
386 3322 : llm = llm + 1
387 : END DO
388 : END DO
389 : END DO
390 : CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
391 618 : ngrd, llm, spdim)
392 :
393 : ! repulsive potential
394 618 : dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
395 618 : dftb_potential(ikind, jkind)%urep_cut = uwork(10)
396 6798 : dftb_potential(ikind, jkind)%urep(:) = 0._dp
397 618 : dftb_potential(ikind, jkind)%urep(1) = uwork(10)
398 1526 : dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
399 :
400 : ! Slater-Koster tables
401 618 : dftb_potential(ikind, jkind)%dgrd = dgrd
402 618 : CALL skreorder(fmat, lmax_a, lmax_b)
403 923578 : dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
404 618 : CALL skreorder(smat, lmax_a, lmax_b)
405 923578 : dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
406 618 : dftb_potential(ikind, jkind)%ngrdcut = ngrd + INT(slako_d0/dgrd)
407 : ! Splines
408 618 : IF (spdim > 0) THEN
409 476 : dftb_potential(ikind, jkind)%s_cut = s_cut
410 1904 : dftb_potential(ikind, jkind)%srep = srep
411 32316 : dftb_potential(ikind, jkind)%spxr = spxr
412 64156 : dftb_potential(ikind, jkind)%scoeff = scoeff
413 1428 : dftb_potential(ikind, jkind)%surr = surr
414 : END IF
415 :
416 618 : DEALLOCATE (fmat)
417 618 : DEALLOCATE (smat)
418 2764 : IF (spdim > 0) THEN
419 476 : DEALLOCATE (spxr)
420 476 : DEALLOCATE (scoeff)
421 : END IF
422 :
423 : END DO
424 :
425 : ! no all other pairs
426 910 : DO ikind = 1, nkind
427 618 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
428 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
429 :
430 618 : IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
431 0 : CALL allocate_dftb_atom_param(dftb_atom_a)
432 0 : CALL set_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
433 : END IF
434 :
435 2950 : DO jkind = 1, nkind
436 :
437 1422 : IF (ikind == jkind) CYCLE
438 804 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
439 804 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
440 :
441 804 : IF (output_unit > 0) THEN
442 0 : WRITE (output_unit, "(A,T30,A50)") " DFTB| Reading parameter file ", &
443 0 : ADJUSTR(TRIM(sk_files(ikind, jkind)))
444 : END IF
445 804 : skself = 0._dp
446 804 : eta = 0._dp
447 804 : occupation = 0._dp
448 804 : IF (para_env%is_source()) THEN
449 402 : runit = get_unit_number()
450 402 : CALL open_file(file_name=sk_files(ikind, jkind), unit_number=runit)
451 : ! grid density and number of grid poin ts
452 402 : 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 402 : 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 402 : READ (runit, fmt=*, END=1, err=1) uwork(1:10)
466 402 : n_urpoly = 0
467 4020 : IF (DOT_PRODUCT(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
468 86 : n_urpoly = 1
469 774 : DO k = 2, 9
470 774 : 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 402 : IF (n_urpoly < 2) n_urpoly = 0
476 : END IF
477 :
478 804 : CALL para_env%bcast(n_urpoly)
479 804 : CALL para_env%bcast(uwork)
480 804 : CALL para_env%bcast(ngrd)
481 804 : CALL para_env%bcast(dgrd)
482 :
483 : ! Slater-Koster table
484 2412 : ALLOCATE (fmat(ngrd, 10))
485 1608 : ALLOCATE (smat(ngrd, 10))
486 804 : IF (para_env%is_source()) THEN
487 199960 : DO k = 1, ngrd
488 199558 : READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
489 2195138 : fmat(k, 1:10) = fwork(1:10)
490 2195540 : smat(k, 1:10) = swork(1:10)
491 : END DO
492 : END IF
493 804 : CALL para_env%bcast(fmat)
494 804 : CALL para_env%bcast(smat)
495 :
496 804 : spdim = 0
497 804 : IF (n_urpoly == 0) THEN
498 632 : IF (para_env%is_source()) THEN
499 : ! Look for spline representation of repulsive potential
500 : search = .TRUE.
501 : DO WHILE (search)
502 6444 : READ (runit, fmt='(A6)', END=1, err=1) cspline
503 6444 : IF (cspline == 'Spline') THEN
504 316 : search = .FALSE.
505 : ! spline dimension and left-hand cutoff
506 316 : READ (runit, fmt=*, END=1, err=1) spdim, s_cut
507 948 : ALLOCATE (spxr(spdim, 2))
508 948 : ALLOCATE (scoeff(spdim, 4))
509 : ! e-functions describing left-hand extrapolation
510 316 : READ (runit, fmt=*, END=1, err=1) srep(1:3)
511 10156 : DO isp = 1, spdim - 1
512 : ! location and coefficients of 'normal' spline range
513 10156 : 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 316 : 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 804 : IF (para_env%is_source()) THEN
523 402 : CALL close_file(unit_number=runit)
524 : END IF
525 :
526 804 : CALL para_env%bcast(spdim)
527 804 : IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
528 948 : ALLOCATE (spxr(spdim, 2))
529 948 : ALLOCATE (scoeff(spdim, 4))
530 : END IF
531 804 : IF (spdim > 0) THEN
532 632 : CALL para_env%bcast(spxr)
533 632 : CALL para_env%bcast(scoeff)
534 632 : CALL para_env%bcast(surr)
535 632 : CALL para_env%bcast(srep)
536 632 : CALL para_env%bcast(s_cut)
537 : END IF
538 :
539 : ! store potential data
540 : ! allocate data
541 804 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
542 804 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
543 804 : llm = 0
544 2428 : DO l1 = 0, MAX(lmax_a, lmax_b)
545 4260 : DO l2 = 0, MIN(l1, lmax_a, lmax_b)
546 5504 : DO m = 0, l2
547 3880 : llm = llm + 1
548 : END DO
549 : END DO
550 : END DO
551 : CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
552 804 : ngrd, llm, spdim)
553 :
554 : ! repulsive potential
555 804 : dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
556 804 : dftb_potential(ikind, jkind)%urep_cut = uwork(10)
557 8844 : dftb_potential(ikind, jkind)%urep(:) = 0._dp
558 804 : dftb_potential(ikind, jkind)%urep(1) = uwork(10)
559 1748 : dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
560 :
561 : ! Slater-Koster tables
562 804 : dftb_potential(ikind, jkind)%dgrd = dgrd
563 804 : CALL skreorder(fmat, lmax_a, lmax_b)
564 1012644 : dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
565 804 : CALL skreorder(smat, lmax_a, lmax_b)
566 1012644 : dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
567 804 : dftb_potential(ikind, jkind)%ngrdcut = ngrd + INT(slako_d0/dgrd)
568 : ! Splines
569 804 : IF (spdim > 0) THEN
570 632 : dftb_potential(ikind, jkind)%s_cut = s_cut
571 2528 : dftb_potential(ikind, jkind)%srep = srep
572 42520 : dftb_potential(ikind, jkind)%spxr = spxr
573 84408 : dftb_potential(ikind, jkind)%scoeff = scoeff
574 1896 : dftb_potential(ikind, jkind)%surr = surr
575 : END IF
576 :
577 804 : DEALLOCATE (fmat)
578 804 : DEALLOCATE (smat)
579 2226 : IF (spdim > 0) THEN
580 632 : DEALLOCATE (spxr)
581 632 : DEALLOCATE (scoeff)
582 : END IF
583 :
584 : END DO
585 : END DO
586 :
587 292 : DEALLOCATE (sk_files)
588 :
589 : ! read dispersion parameters (UFF type)
590 292 : IF (dftb_control%dispersion) THEN
591 :
592 112 : IF (dftb_control%dispersion_type == dispersion_uff) THEN
593 : file_name = TRIM(dftb_control%sk_file_path)//"/"// &
594 94 : TRIM(dftb_control%uff_force_field)
595 : BLOCK
596 : TYPE(cp_parser_type) :: parser
597 476 : DO ikind = 1, nkind
598 194 : CALL get_atomic_kind(atomic_kind_set(ikind), name=iname)
599 194 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
600 :
601 194 : m = LEN_TRIM(iname)
602 194 : CALL parser_create(parser, file_name, para_env=para_env)
603 194 : found = .FALSE.
604 : DO
605 : at_end = .FALSE.
606 2104 : CALL parser_get_next_line(parser, 1, at_end)
607 2104 : IF (at_end) EXIT
608 2104 : CALL parser_get_object(parser, name_a)
609 : ! parser is no longer removing leading quotes
610 2104 : IF (name_a(1:1) == '"') name_a(1:m) = name_a(2:m + 1)
611 2104 : IF (name_a(1:m) == TRIM(iname)) THEN
612 194 : CALL parser_get_object(parser, rb)
613 194 : CALL parser_get_object(parser, rb)
614 194 : CALL parser_get_object(parser, ra)
615 194 : CALL parser_get_object(parser, da)
616 194 : found = .TRUE.
617 194 : ra = ra/angstrom
618 194 : da = da/kcalmol
619 194 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, name=iname, xi=ra, di=da)
620 194 : EXIT
621 : END IF
622 : END DO
623 482 : 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 910 : DO ikind = 1, nkind
632 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
633 : radmax = (dftb_potential(ikind, ikind)%ngrdcut + 1)* &
634 618 : dftb_potential(ikind, ikind)%dgrd*0.5_dp
635 910 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=radmax)
636 : END DO
637 910 : DO ikind = 1, nkind
638 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
639 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
640 2332 : DO jkind = 1, nkind
641 1422 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
642 1422 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, cutoff=rb)
643 : radmax = (dftb_potential(ikind, jkind)%ngrdcut + 1)* &
644 1422 : dftb_potential(ikind, jkind)%dgrd
645 2040 : IF (ra + rb < radmax) THEN
646 8 : ra = ra + (radmax - ra - rb)*0.5_dp
647 8 : rb = rb + (radmax - ra - rb)*0.5_dp
648 8 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
649 8 : 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 910 : DO ikind = 1, nkind
656 618 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
657 618 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, zeff=zeff)
658 : CALL set_potential(potential=qs_kind_set(ikind)%all_potential, &
659 910 : zeff=zeff, zeff_correction=0.0_dp)
660 : END DO
661 :
662 : ! setup DFTB3 parameters
663 292 : IF (dftb_control%dftb3_diagonal) THEN
664 134 : DO ikind = 1, nkind
665 88 : CALL get_qs_kind(qs_kind_set(ikind), dftb3_param=db)
666 88 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
667 222 : 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 292 : IF (dftb_control%dispersion) THEN
673 112 : IF (dftb_control%dispersion_type == dispersion_uff) THEN
674 94 : eps_disp = dftb_control%eps_disp
675 288 : DO ikind = 1, nkind
676 194 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
677 194 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, xi=ra, di=da)
678 194 : rcdisp = 0._dp
679 652 : DO jkind = 1, nkind
680 458 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
681 458 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, xi=rb, di=db)
682 458 : xij = SQRT(ra*rb)
683 458 : dij = SQRT(da*db)
684 458 : dftb_potential(ikind, jkind)%xij = xij
685 458 : dftb_potential(ikind, jkind)%dij = dij
686 458 : dftb_potential(ikind, jkind)%x0ij = xij*(0.5_dp**(1.0_dp/6.0_dp))
687 458 : dftb_potential(ikind, jkind)%a = dij*396.0_dp/25.0_dp
688 : dftb_potential(ikind, jkind)%b = &
689 458 : 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 458 : -dij/(xij**10)*2.0_dp**(2.0_dp/3.0_dp)*552.0_dp/25.0_dp
692 458 : rmax6 = ((8._dp*pi*dij/eps_disp)*xij**6)**0.25_dp
693 652 : rcdisp = MAX(rcdisp, rmax6*0.5_dp)
694 : END DO
695 288 : 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 :
702 : 1 CONTINUE
703 0 : CPABORT("")
704 :
705 292 : 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 2844 : 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 1411440 : DO i = 1, SIZE(xmat, 1)
738 1408596 : skllm = 0._dp
739 1408596 : skllm(0, 0, 0) = xmat(i, 10)
740 1408596 : skllm(1, 0, 0) = xmat(i, 9)
741 1408596 : skllm(2, 0, 0) = xmat(i, 8)
742 1408596 : skllm(1, 1, 1) = xmat(i, 7)
743 1408596 : skllm(1, 1, 0) = xmat(i, 6)
744 1408596 : skllm(2, 1, 1) = xmat(i, 5)
745 1408596 : skllm(2, 1, 0) = xmat(i, 4)
746 1408596 : skllm(2, 2, 2) = xmat(i, 3)
747 1408596 : skllm(2, 2, 1) = xmat(i, 2)
748 1408596 : skllm(2, 2, 0) = xmat(i, 1)
749 1408596 : llm = 0
750 4018800 : DO l1 = 0, MAX(la, lb)
751 7234776 : DO l2 = 0, MIN(l1, la, lb)
752 9687916 : DO m = 0, l2
753 3861736 : llm = llm + 1
754 7080556 : xmat(i, llm) = skllm(l1, l2, m)
755 : END DO
756 : END DO
757 : END DO
758 : END DO
759 : !
760 2844 : END SUBROUTINE skreorder
761 :
762 : END MODULE qs_dftb_parameters
763 :
|
