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