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 254 : 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 254 : 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 254 : 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 254 : output_unit = -1
112 254 : NULLIFY (logger)
113 254 : logger => cp_get_default_logger()
114 254 : 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 254 : sklist = (dftb_control%sk_file_list /= "")
145 :
146 254 : nkind = SIZE(atomic_kind_set)
147 1016 : ALLOCATE (sk_files(nkind, nkind))
148 : ! allocate potential structures
149 6582 : ALLOCATE (dftb_potential(nkind, nkind))
150 254 : CALL qs_dftb_pairpot_init(dftb_potential)
151 :
152 788 : DO ikind = 1, nkind
153 534 : CALL get_atomic_kind(atomic_kind_set(ikind), name=iname, element_symbol=iel)
154 534 : CALL uppercase(iname)
155 534 : CALL uppercase(iel)
156 534 : ldum = qmmm_ff_precond_only_qm(iname)
157 2010 : DO jkind = 1, nkind
158 1222 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname, element_symbol=jel)
159 1222 : CALL uppercase(jname)
160 1222 : CALL uppercase(jel)
161 1222 : ldum = qmmm_ff_precond_only_qm(jname)
162 1222 : found = .FALSE.
163 1306 : DO k = 1, SIZE(dftb_control%sk_pair_list, 2)
164 170 : name_a = TRIM(dftb_control%sk_pair_list(1, k))
165 170 : name_b = TRIM(dftb_control%sk_pair_list(2, k))
166 170 : CALL uppercase(name_a)
167 170 : CALL uppercase(name_b)
168 1306 : IF ((iname == name_a .AND. jname == name_b)) THEN
169 : sk_files(ikind, jkind) = TRIM(dftb_control%sk_file_path)//"/"// &
170 86 : TRIM(dftb_control%sk_pair_list(3, k))
171 86 : found = .TRUE.
172 86 : EXIT
173 : END IF
174 : END DO
175 1222 : IF (.NOT. found .AND. sklist) THEN
176 : file_name = TRIM(dftb_control%sk_file_path)//"/"// &
177 1136 : TRIM(dftb_control%sk_file_list)
178 1136 : BLOCK
179 : TYPE(cp_parser_type) :: parser
180 1136 : CALL parser_create(parser, file_name, para_env=para_env)
181 : DO
182 : at_end = .FALSE.
183 16832 : CALL parser_get_next_line(parser, 1, at_end)
184 16832 : IF (at_end) EXIT
185 16832 : CALL parser_get_object(parser, name_a, lower_to_upper=.TRUE.)
186 16832 : CALL parser_get_object(parser, name_b, lower_to_upper=.TRUE.)
187 : !Checking Names
188 16832 : IF ((iname == name_a .AND. jname == name_b)) THEN
189 1136 : CALL parser_get_object(parser, skfn, string_length=8)
190 : sk_files(ikind, jkind) = TRIM(dftb_control%sk_file_path)//"/"// &
191 1136 : TRIM(skfn)
192 1136 : found = .TRUE.
193 1136 : EXIT
194 : END IF
195 : !Checking Element
196 15696 : 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 4544 : CALL parser_release(parser)
205 : END BLOCK
206 : END IF
207 1222 : IF (.NOT. found) &
208 : CALL cp_abort(__LOCATION__, &
209 : "Failure in assigning KINDS <"//TRIM(iname)//"> and <"//TRIM(jname)// &
210 534 : "> to a DFTB interaction pair!")
211 : END DO
212 : END DO
213 : ! reading the files
214 : ! read all pairs, equal kind first
215 788 : DO ikind = 1, nkind
216 534 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
217 :
218 534 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
219 534 : IF (.NOT. ASSOCIATED(dftb_atom_a)) THEN
220 534 : CALL allocate_dftb_atom_param(dftb_atom_a)
221 534 : 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 534 : jkind = ikind
226 :
227 534 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
228 534 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
229 :
230 534 : 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 534 : skself = 0._dp
235 534 : eta = 0._dp
236 534 : occupation = 0._dp
237 534 : IF (para_env%is_source()) THEN
238 267 : runit = get_unit_number()
239 267 : CALL open_file(file_name=sk_files(jkind, ikind), unit_number=runit)
240 : ! grid density and number of grid poin ts
241 267 : 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 267 : ngrd = ngrd - 1
247 : !
248 : ! orbital energy, total energy, hardness, occupation
249 267 : READ (runit, fmt=*, END=1, err=1) skself(2:0:-1), energy, &
250 534 : eta(2:0:-1), occupation(2:0:-1)
251 : ! repulsive potential as polynomial
252 267 : READ (runit, fmt=*, END=1, err=1) uwork(1:10)
253 267 : n_urpoly = 0
254 2670 : IF (DOT_PRODUCT(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
255 62 : n_urpoly = 1
256 558 : DO k = 2, 9
257 558 : 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 267 : IF (n_urpoly < 2) n_urpoly = 0
263 : END IF
264 :
265 534 : CALL para_env%bcast(n_urpoly)
266 534 : CALL para_env%bcast(uwork)
267 534 : CALL para_env%bcast(ngrd)
268 534 : CALL para_env%bcast(dgrd)
269 :
270 534 : CALL para_env%bcast(skself)
271 534 : CALL para_env%bcast(energy)
272 534 : CALL para_env%bcast(eta)
273 534 : 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 2670 : skself=skself, energy=energy, eta=eta, occupation=occupation)
278 :
279 : ! Slater-Koster table
280 1602 : ALLOCATE (fmat(ngrd, 10))
281 1068 : ALLOCATE (smat(ngrd, 10))
282 534 : IF (para_env%is_source()) THEN
283 131800 : DO k = 1, ngrd
284 131533 : READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
285 1446863 : fmat(k, 1:10) = fwork(1:10)
286 1447130 : smat(k, 1:10) = swork(1:10)
287 : END DO
288 : END IF
289 534 : CALL para_env%bcast(fmat)
290 534 : 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 534 : CALL get_qs_kind(qs_kind_set(ikind), lmax_dftb=lmax)
304 534 : IF (lmax < 0) THEN
305 530 : lmax = 0
306 2650 : DO l = 0, 3
307 0 : SELECT CASE (l)
308 : CASE DEFAULT
309 0 : CPABORT("")
310 : CASE (0)
311 530 : lp = 10
312 : CASE (1)
313 530 : lp = 7
314 : CASE (2)
315 530 : lp = 3
316 : CASE (3)
317 2120 : 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 946000 : IF ((ABS(skself(l)) > 0._dp) .OR. &
322 1418 : (SUM(ABS(fmat(ngrd/10:ngrd, lp))) > 0._dp)) lmax = l
323 : END DO
324 : ! l=2 (d) is maximum
325 530 : lmax = MIN(2, lmax)
326 : END IF
327 534 : 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 534 : lmax=lmax, natorb=(lmax + 1)**2)
334 :
335 534 : spdim = 0
336 534 : IF (n_urpoly == 0) THEN
337 410 : IF (para_env%is_source()) THEN
338 : ! Look for spline representation of repulsive potential
339 : search = .TRUE.
340 : DO WHILE (search)
341 4306 : READ (runit, fmt='(A6)', END=1, err=1) cspline
342 4306 : IF (cspline == 'Spline') THEN
343 205 : search = .FALSE.
344 : ! spline dimension and left-hand cutoff
345 205 : READ (runit, fmt=*, END=1, err=1) spdim, s_cut
346 615 : ALLOCATE (spxr(spdim, 2))
347 615 : ALLOCATE (scoeff(spdim, 4))
348 : ! e-functions describing left-hand extrapolation
349 205 : READ (runit, fmt=*, END=1, err=1) srep(1:3)
350 6627 : DO isp = 1, spdim - 1
351 : ! location and coefficients of 'normal' spline range
352 6627 : 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 205 : 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 534 : IF (para_env%is_source()) THEN
362 267 : CALL close_file(unit_number=runit)
363 : END IF
364 :
365 534 : CALL para_env%bcast(spdim)
366 534 : IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
367 615 : ALLOCATE (spxr(spdim, 2))
368 615 : ALLOCATE (scoeff(spdim, 4))
369 : END IF
370 534 : IF (spdim > 0) THEN
371 410 : CALL para_env%bcast(spxr)
372 410 : CALL para_env%bcast(scoeff)
373 410 : CALL para_env%bcast(surr)
374 410 : CALL para_env%bcast(srep)
375 410 : CALL para_env%bcast(s_cut)
376 : END IF
377 :
378 : ! store potential data
379 : ! allocate data
380 534 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
381 534 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
382 534 : llm = 0
383 1410 : DO l1 = 0, MAX(lmax_a, lmax_b)
384 2648 : DO l2 = 0, MIN(l1, lmax_a, lmax_b)
385 3734 : DO m = 0, l2
386 2858 : llm = llm + 1
387 : END DO
388 : END DO
389 : END DO
390 : CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
391 534 : ngrd, llm, spdim)
392 :
393 : ! repulsive potential
394 534 : dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
395 534 : dftb_potential(ikind, jkind)%urep_cut = uwork(10)
396 5874 : dftb_potential(ikind, jkind)%urep(:) = 0._dp
397 534 : dftb_potential(ikind, jkind)%urep(1) = uwork(10)
398 1328 : dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
399 :
400 : ! Slater-Koster tables
401 534 : dftb_potential(ikind, jkind)%dgrd = dgrd
402 534 : CALL skreorder(fmat, lmax_a, lmax_b)
403 789494 : dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
404 534 : CALL skreorder(smat, lmax_a, lmax_b)
405 789494 : dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
406 534 : dftb_potential(ikind, jkind)%ngrdcut = ngrd + INT(slako_d0/dgrd)
407 : ! Splines
408 534 : IF (spdim > 0) THEN
409 410 : dftb_potential(ikind, jkind)%s_cut = s_cut
410 1640 : dftb_potential(ikind, jkind)%srep = srep
411 27738 : dftb_potential(ikind, jkind)%spxr = spxr
412 55066 : dftb_potential(ikind, jkind)%scoeff = scoeff
413 1230 : dftb_potential(ikind, jkind)%surr = surr
414 : END IF
415 :
416 534 : DEALLOCATE (fmat)
417 534 : DEALLOCATE (smat)
418 2390 : IF (spdim > 0) THEN
419 410 : DEALLOCATE (spxr)
420 410 : DEALLOCATE (scoeff)
421 : END IF
422 :
423 : END DO
424 :
425 : ! no all other pairs
426 788 : DO ikind = 1, nkind
427 534 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=iname)
428 534 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
429 :
430 534 : 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 2544 : DO jkind = 1, nkind
436 :
437 1222 : IF (ikind == jkind) CYCLE
438 688 : CALL get_atomic_kind(atomic_kind_set(jkind), name=jname)
439 688 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
440 :
441 688 : 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 688 : skself = 0._dp
446 688 : eta = 0._dp
447 688 : occupation = 0._dp
448 688 : IF (para_env%is_source()) THEN
449 344 : runit = get_unit_number()
450 344 : CALL open_file(file_name=sk_files(ikind, jkind), unit_number=runit)
451 : ! grid density and number of grid poin ts
452 344 : 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 344 : 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 344 : READ (runit, fmt=*, END=1, err=1) uwork(1:10)
466 344 : n_urpoly = 0
467 3440 : IF (DOT_PRODUCT(uwork(2:10), uwork(2:10)) >= 1.e-12_dp) THEN
468 68 : n_urpoly = 1
469 612 : DO k = 2, 9
470 612 : 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 344 : IF (n_urpoly < 2) n_urpoly = 0
476 : END IF
477 :
478 688 : CALL para_env%bcast(n_urpoly)
479 688 : CALL para_env%bcast(uwork)
480 688 : CALL para_env%bcast(ngrd)
481 688 : CALL para_env%bcast(dgrd)
482 :
483 : ! Slater-Koster table
484 2064 : ALLOCATE (fmat(ngrd, 10))
485 1376 : ALLOCATE (smat(ngrd, 10))
486 688 : IF (para_env%is_source()) THEN
487 170760 : DO k = 1, ngrd
488 170416 : READ (runit, fmt=*, END=1, err=1) fwork(1:10), swork(1:10)
489 1874576 : fmat(k, 1:10) = fwork(1:10)
490 1874920 : smat(k, 1:10) = swork(1:10)
491 : END DO
492 : END IF
493 688 : CALL para_env%bcast(fmat)
494 688 : CALL para_env%bcast(smat)
495 :
496 688 : spdim = 0
497 688 : IF (n_urpoly == 0) THEN
498 552 : IF (para_env%is_source()) THEN
499 : ! Look for spline representation of repulsive potential
500 : search = .TRUE.
501 : DO WHILE (search)
502 5630 : READ (runit, fmt='(A6)', END=1, err=1) cspline
503 5630 : IF (cspline == 'Spline') THEN
504 276 : search = .FALSE.
505 : ! spline dimension and left-hand cutoff
506 276 : READ (runit, fmt=*, END=1, err=1) spdim, s_cut
507 828 : ALLOCATE (spxr(spdim, 2))
508 828 : ALLOCATE (scoeff(spdim, 4))
509 : ! e-functions describing left-hand extrapolation
510 276 : READ (runit, fmt=*, END=1, err=1) srep(1:3)
511 8772 : DO isp = 1, spdim - 1
512 : ! location and coefficients of 'normal' spline range
513 8772 : 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 276 : 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 688 : IF (para_env%is_source()) THEN
523 344 : CALL close_file(unit_number=runit)
524 : END IF
525 :
526 688 : CALL para_env%bcast(spdim)
527 688 : IF (spdim > 0 .AND. (.NOT. para_env%is_source())) THEN
528 828 : ALLOCATE (spxr(spdim, 2))
529 828 : ALLOCATE (scoeff(spdim, 4))
530 : END IF
531 688 : IF (spdim > 0) THEN
532 552 : CALL para_env%bcast(spxr)
533 552 : CALL para_env%bcast(scoeff)
534 552 : CALL para_env%bcast(surr)
535 552 : CALL para_env%bcast(srep)
536 552 : CALL para_env%bcast(s_cut)
537 : END IF
538 :
539 : ! store potential data
540 : ! allocate data
541 688 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, lmax=lmax_a)
542 688 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, lmax=lmax_b)
543 688 : llm = 0
544 2076 : DO l1 = 0, MAX(lmax_a, lmax_b)
545 3636 : DO l2 = 0, MIN(l1, lmax_a, lmax_b)
546 4684 : DO m = 0, l2
547 3296 : llm = llm + 1
548 : END DO
549 : END DO
550 : END DO
551 : CALL qs_dftb_pairpot_create(dftb_potential(ikind, jkind), &
552 688 : ngrd, llm, spdim)
553 :
554 : ! repulsive potential
555 688 : dftb_potential(ikind, jkind)%n_urpoly = n_urpoly
556 688 : dftb_potential(ikind, jkind)%urep_cut = uwork(10)
557 7568 : dftb_potential(ikind, jkind)%urep(:) = 0._dp
558 688 : dftb_potential(ikind, jkind)%urep(1) = uwork(10)
559 1440 : dftb_potential(ikind, jkind)%urep(2:n_urpoly) = uwork(2:n_urpoly)
560 :
561 : ! Slater-Koster tables
562 688 : dftb_potential(ikind, jkind)%dgrd = dgrd
563 688 : CALL skreorder(fmat, lmax_a, lmax_b)
564 852528 : dftb_potential(ikind, jkind)%fmat(:, 1:llm) = fmat(:, 1:llm)
565 688 : CALL skreorder(smat, lmax_a, lmax_b)
566 852528 : dftb_potential(ikind, jkind)%smat(:, 1:llm) = smat(:, 1:llm)
567 688 : dftb_potential(ikind, jkind)%ngrdcut = ngrd + INT(slako_d0/dgrd)
568 : ! Splines
569 688 : IF (spdim > 0) THEN
570 552 : dftb_potential(ikind, jkind)%s_cut = s_cut
571 2208 : dftb_potential(ikind, jkind)%srep = srep
572 36744 : dftb_potential(ikind, jkind)%spxr = spxr
573 72936 : dftb_potential(ikind, jkind)%scoeff = scoeff
574 1656 : dftb_potential(ikind, jkind)%surr = surr
575 : END IF
576 :
577 688 : DEALLOCATE (fmat)
578 688 : DEALLOCATE (smat)
579 1910 : IF (spdim > 0) THEN
580 552 : DEALLOCATE (spxr)
581 552 : DEALLOCATE (scoeff)
582 : END IF
583 :
584 : END DO
585 : END DO
586 :
587 254 : DEALLOCATE (sk_files)
588 :
589 : ! read dispersion parameters (UFF type)
590 254 : IF (dftb_control%dispersion) THEN
591 :
592 106 : IF (dftb_control%dispersion_type == dispersion_uff) THEN
593 : file_name = TRIM(dftb_control%sk_file_path)//"/"// &
594 88 : TRIM(dftb_control%uff_force_field)
595 : BLOCK
596 : TYPE(cp_parser_type) :: parser
597 444 : DO ikind = 1, nkind
598 180 : CALL get_atomic_kind(atomic_kind_set(ikind), name=iname)
599 180 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
600 :
601 180 : m = LEN_TRIM(iname)
602 180 : CALL parser_create(parser, file_name, para_env=para_env)
603 180 : found = .FALSE.
604 : DO
605 : at_end = .FALSE.
606 1988 : CALL parser_get_next_line(parser, 1, at_end)
607 1988 : IF (at_end) EXIT
608 1988 : CALL parser_get_object(parser, name_a)
609 : ! parser is no longer removing leading quotes
610 1988 : IF (name_a(1:1) == '"') name_a(1:m) = name_a(2:m + 1)
611 1988 : IF (name_a(1:m) == TRIM(iname)) THEN
612 180 : CALL parser_get_object(parser, rb)
613 180 : CALL parser_get_object(parser, rb)
614 180 : CALL parser_get_object(parser, ra)
615 180 : CALL parser_get_object(parser, da)
616 180 : found = .TRUE.
617 180 : ra = ra/angstrom
618 180 : da = da/kcalmol
619 180 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, name=iname, xi=ra, di=da)
620 180 : EXIT
621 : END IF
622 : END DO
623 448 : 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 788 : DO ikind = 1, nkind
632 534 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
633 : radmax = (dftb_potential(ikind, ikind)%ngrdcut + 1)* &
634 534 : dftb_potential(ikind, ikind)%dgrd*0.5_dp
635 788 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=radmax)
636 : END DO
637 788 : DO ikind = 1, nkind
638 534 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
639 534 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
640 2010 : DO jkind = 1, nkind
641 1222 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
642 1222 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, cutoff=rb)
643 : radmax = (dftb_potential(ikind, jkind)%ngrdcut + 1)* &
644 1222 : dftb_potential(ikind, jkind)%dgrd
645 1756 : IF (ra + rb < radmax) THEN
646 4 : ra = ra + (radmax - ra - rb)*0.5_dp
647 4 : rb = rb + (radmax - ra - rb)*0.5_dp
648 4 : CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a, cutoff=ra)
649 4 : 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 788 : DO ikind = 1, nkind
656 534 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
657 534 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, zeff=zeff)
658 : CALL set_potential(potential=qs_kind_set(ikind)%all_potential, &
659 788 : zeff=zeff, zeff_correction=0.0_dp)
660 : END DO
661 :
662 : ! setup DFTB3 parameters
663 254 : IF (dftb_control%dftb3_diagonal) THEN
664 88 : DO ikind = 1, nkind
665 58 : CALL get_qs_kind(qs_kind_set(ikind), dftb3_param=db)
666 58 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
667 146 : 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 254 : IF (dftb_control%dispersion) THEN
673 106 : IF (dftb_control%dispersion_type == dispersion_uff) THEN
674 88 : eps_disp = dftb_control%eps_disp
675 268 : DO ikind = 1, nkind
676 180 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom_a)
677 180 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_a, xi=ra, di=da)
678 180 : rcdisp = 0._dp
679 600 : DO jkind = 1, nkind
680 420 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_atom_b)
681 420 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom_b, xi=rb, di=db)
682 420 : xij = SQRT(ra*rb)
683 420 : dij = SQRT(da*db)
684 420 : dftb_potential(ikind, jkind)%xij = xij
685 420 : dftb_potential(ikind, jkind)%dij = dij
686 420 : dftb_potential(ikind, jkind)%x0ij = xij*(0.5_dp**(1.0_dp/6.0_dp))
687 420 : dftb_potential(ikind, jkind)%a = dij*396.0_dp/25.0_dp
688 : dftb_potential(ikind, jkind)%b = &
689 420 : 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 420 : -dij/(xij**10)*2.0_dp**(2.0_dp/3.0_dp)*552.0_dp/25.0_dp
692 420 : rmax6 = ((8._dp*pi*dij/eps_disp)*xij**6)**0.25_dp
693 600 : rcdisp = MAX(rcdisp, rmax6*0.5_dp)
694 : END DO
695 268 : 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 254 : 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 2444 : 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 1210240 : DO i = 1, SIZE(xmat, 1)
738 1207796 : skllm = 0._dp
739 1207796 : skllm(0, 0, 0) = xmat(i, 10)
740 1207796 : skllm(1, 0, 0) = xmat(i, 9)
741 1207796 : skllm(2, 0, 0) = xmat(i, 8)
742 1207796 : skllm(1, 1, 1) = xmat(i, 7)
743 1207796 : skllm(1, 1, 0) = xmat(i, 6)
744 1207796 : skllm(2, 1, 1) = xmat(i, 5)
745 1207796 : skllm(2, 1, 0) = xmat(i, 4)
746 1207796 : skllm(2, 2, 2) = xmat(i, 3)
747 1207796 : skllm(2, 2, 1) = xmat(i, 2)
748 1207796 : skllm(2, 2, 0) = xmat(i, 1)
749 1207796 : llm = 0
750 3438752 : DO l1 = 0, MAX(la, lb)
751 6176872 : DO l2 = 0, MIN(l1, la, lb)
752 8243964 : DO m = 0, l2
753 3274888 : llm = llm + 1
754 6015452 : xmat(i, llm) = skllm(l1, l2, m)
755 : END DO
756 : END DO
757 : END DO
758 : END DO
759 : !
760 2444 : END SUBROUTINE skreorder
761 :
762 : END MODULE qs_dftb_parameters
763 :
|
