(git:d3d49ac)
Loading...
Searching...
No Matches
qs_dftb_utils.F
Go to the documentation of this file.
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!> \brief Working with the DFTB parameter types.
10!> \author JGH (24.02.2007)
11! **************************************************************************************************
12MODULE qs_dftb_utils
13
14 USE cp_log_handling, ONLY: cp_get_default_logger,&
15 cp_logger_type
16 USE cp_output_handling, ONLY: cp_p_file,&
17 cp_print_key_finished_output,&
18 cp_print_key_should_output,&
19 cp_print_key_unit_nr
20 USE input_section_types, ONLY: section_vals_type
21 USE kinds, ONLY: default_string_length,&
22 dp
23 USE qs_dftb_types, ONLY: qs_dftb_atom_type
24#include "./base/base_uses.f90"
25
26 IMPLICIT NONE
27
28 PRIVATE
29
30 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_utils'
31
32 ! Maximum number of points used for interpolation
33 INTEGER, PARAMETER :: max_inter = 5
34 ! Maximum number of points used for extrapolation
35 INTEGER, PARAMETER :: max_extra = 9
36 ! see also qs_dftb_parameters
37 REAL(dp), PARAMETER :: slako_d0 = 1._dp
38 ! pointer to skab
39 INTEGER, DIMENSION(0:3, 0:3, 0:3, 0:3, 0:3):: iptr
40 ! small real number
41 REAL(dp), PARAMETER :: rtiny = 1.e-10_dp
42 ! eta(0) for mm atoms and non-scc qm atoms
43 REAL(dp), PARAMETER :: eta_mm = 0.47_dp
44 ! step size for qmmm finite difference
45 REAL(dp), PARAMETER :: ddrmm = 0.0001_dp
46
47 PUBLIC :: allocate_dftb_atom_param, &
48 deallocate_dftb_atom_param, &
49 get_dftb_atom_param, &
50 set_dftb_atom_param, &
51 write_dftb_atom_param
52 PUBLIC :: compute_block_sk, &
53 urep_egr, iptr
54
55CONTAINS
56
57! **************************************************************************************************
58!> \brief ...
59!> \param dftb_parameter ...
60! **************************************************************************************************
61 SUBROUTINE allocate_dftb_atom_param(dftb_parameter)
62
63 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
64
65 IF (ASSOCIATED(dftb_parameter)) THEN
66 CALL deallocate_dftb_atom_param(dftb_parameter)
67 END IF
68
69 ALLOCATE (dftb_parameter)
70
71 dftb_parameter%defined = .false.
72 dftb_parameter%name = ""
73 dftb_parameter%typ = "NONE"
74 dftb_parameter%z = -1
75 dftb_parameter%zeff = -1.0_dp
76 dftb_parameter%natorb = 0
77 dftb_parameter%lmax = -1
78 dftb_parameter%skself = 0.0_dp
79 dftb_parameter%occupation = 0.0_dp
80 dftb_parameter%eta = 0.0_dp
81 dftb_parameter%energy = 0.0_dp
82 dftb_parameter%xi = 0.0_dp
83 dftb_parameter%di = 0.0_dp
84 dftb_parameter%rcdisp = 0.0_dp
85 dftb_parameter%dudq = 0.0_dp
86
87 END SUBROUTINE allocate_dftb_atom_param
88
89! **************************************************************************************************
90!> \brief ...
91!> \param dftb_parameter ...
92! **************************************************************************************************
93 SUBROUTINE deallocate_dftb_atom_param(dftb_parameter)
94
95 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
96
97 cpassert(ASSOCIATED(dftb_parameter))
98 DEALLOCATE (dftb_parameter)
99
100 END SUBROUTINE deallocate_dftb_atom_param
101
102! **************************************************************************************************
103!> \brief ...
104!> \param dftb_parameter ...
105!> \param name ...
106!> \param typ ...
107!> \param defined ...
108!> \param z ...
109!> \param zeff ...
110!> \param natorb ...
111!> \param lmax ...
112!> \param skself ...
113!> \param occupation ...
114!> \param eta ...
115!> \param energy ...
116!> \param cutoff ...
117!> \param xi ...
118!> \param di ...
119!> \param rcdisp ...
120!> \param dudq ...
121! **************************************************************************************************
122 SUBROUTINE get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, &
123 lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
124
125 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
126 CHARACTER(LEN=default_string_length), &
127 INTENT(OUT), OPTIONAL :: name, typ
128 LOGICAL, INTENT(OUT), OPTIONAL :: defined
129 INTEGER, INTENT(OUT), OPTIONAL :: z
130 REAL(kind=dp), INTENT(OUT), OPTIONAL :: zeff
131 INTEGER, INTENT(OUT), OPTIONAL :: natorb, lmax
132 REAL(kind=dp), DIMENSION(0:3), OPTIONAL :: skself, occupation, eta
133 REAL(kind=dp), OPTIONAL :: energy, cutoff, xi, di, rcdisp, dudq
134
135 cpassert(ASSOCIATED(dftb_parameter))
136
137 IF (PRESENT(name)) name = dftb_parameter%name
138 IF (PRESENT(typ)) typ = dftb_parameter%typ
139 IF (PRESENT(defined)) defined = dftb_parameter%defined
140 IF (PRESENT(z)) z = dftb_parameter%z
141 IF (PRESENT(zeff)) zeff = dftb_parameter%zeff
142 IF (PRESENT(natorb)) natorb = dftb_parameter%natorb
143 IF (PRESENT(lmax)) lmax = dftb_parameter%lmax
144 IF (PRESENT(skself)) skself = dftb_parameter%skself
145 IF (PRESENT(eta)) eta = dftb_parameter%eta
146 IF (PRESENT(energy)) energy = dftb_parameter%energy
147 IF (PRESENT(cutoff)) cutoff = dftb_parameter%cutoff
148 IF (PRESENT(occupation)) occupation = dftb_parameter%occupation
149 IF (PRESENT(xi)) xi = dftb_parameter%xi
150 IF (PRESENT(di)) di = dftb_parameter%di
151 IF (PRESENT(rcdisp)) rcdisp = dftb_parameter%rcdisp
152 IF (PRESENT(dudq)) dudq = dftb_parameter%dudq
153
154 END SUBROUTINE get_dftb_atom_param
155
156! **************************************************************************************************
157!> \brief ...
158!> \param dftb_parameter ...
159!> \param name ...
160!> \param typ ...
161!> \param defined ...
162!> \param z ...
163!> \param zeff ...
164!> \param natorb ...
165!> \param lmax ...
166!> \param skself ...
167!> \param occupation ...
168!> \param eta ...
169!> \param energy ...
170!> \param cutoff ...
171!> \param xi ...
172!> \param di ...
173!> \param rcdisp ...
174!> \param dudq ...
175! **************************************************************************************************
176 SUBROUTINE set_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, &
177 lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
178
179 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
180 CHARACTER(LEN=default_string_length), INTENT(IN), &
181 OPTIONAL :: name, typ
182 LOGICAL, INTENT(IN), OPTIONAL :: defined
183 INTEGER, INTENT(IN), OPTIONAL :: z
184 REAL(kind=dp), INTENT(IN), OPTIONAL :: zeff
185 INTEGER, INTENT(IN), OPTIONAL :: natorb, lmax
186 REAL(kind=dp), DIMENSION(0:3), OPTIONAL :: skself, occupation, eta
187 REAL(kind=dp), OPTIONAL :: energy, cutoff, xi, di, rcdisp, dudq
188
189 cpassert(ASSOCIATED(dftb_parameter))
190
191 IF (PRESENT(name)) dftb_parameter%name = name
192 IF (PRESENT(typ)) dftb_parameter%typ = typ
193 IF (PRESENT(defined)) dftb_parameter%defined = defined
194 IF (PRESENT(z)) dftb_parameter%z = z
195 IF (PRESENT(zeff)) dftb_parameter%zeff = zeff
196 IF (PRESENT(natorb)) dftb_parameter%natorb = natorb
197 IF (PRESENT(lmax)) dftb_parameter%lmax = lmax
198 IF (PRESENT(skself)) dftb_parameter%skself = skself
199 IF (PRESENT(eta)) dftb_parameter%eta = eta
200 IF (PRESENT(occupation)) dftb_parameter%occupation = occupation
201 IF (PRESENT(energy)) dftb_parameter%energy = energy
202 IF (PRESENT(cutoff)) dftb_parameter%cutoff = cutoff
203 IF (PRESENT(xi)) dftb_parameter%xi = xi
204 IF (PRESENT(di)) dftb_parameter%di = di
205 IF (PRESENT(rcdisp)) dftb_parameter%rcdisp = rcdisp
206 IF (PRESENT(dudq)) dftb_parameter%dudq = dudq
207
208 END SUBROUTINE set_dftb_atom_param
209
210! **************************************************************************************************
211!> \brief ...
212!> \param dftb_parameter ...
213!> \param subsys_section ...
214! **************************************************************************************************
215 SUBROUTINE write_dftb_atom_param(dftb_parameter, subsys_section)
216
217 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
218 TYPE(section_vals_type), POINTER :: subsys_section
219
220 CHARACTER(LEN=default_string_length) :: name, typ
221 INTEGER :: lmax, natorb, output_unit, z
222 LOGICAL :: defined
223 REAL(dp) :: zeff
224 TYPE(cp_logger_type), POINTER :: logger
225
226 NULLIFY (logger)
227 logger => cp_get_default_logger()
228 IF (ASSOCIATED(dftb_parameter) .AND. &
229 btest(cp_print_key_should_output(logger%iter_info, subsys_section, &
230 "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN
231
232 output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", &
233 extension=".Log")
234
235 IF (output_unit > 0) THEN
236 CALL get_dftb_atom_param(dftb_parameter, name=name, typ=typ, defined=defined, &
237 z=z, zeff=zeff, natorb=natorb, lmax=lmax)
238
239 WRITE (unit=output_unit, fmt="(/,A,T67,A14)") &
240 " DFTB parameters: ", trim(name)
241 IF (defined) THEN
242 WRITE (unit=output_unit, fmt="(T16,A,T71,F10.2)") &
243 "Effective core charge:", zeff
244 WRITE (unit=output_unit, fmt="(T16,A,T71,I10)") &
245 "Number of orbitals:", natorb
246 ELSE
247 WRITE (unit=output_unit, fmt="(T55,A)") &
248 "Parameters are not defined"
249 END IF
250 END IF
251 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
252 "PRINT%KINDS")
253 END IF
254
255 END SUBROUTINE write_dftb_atom_param
256
257! **************************************************************************************************
258!> \brief ...
259!> \param block ...
260!> \param smatij ...
261!> \param smatji ...
262!> \param rij ...
263!> \param ngrd ...
264!> \param ngrdcut ...
265!> \param dgrd ...
266!> \param llm ...
267!> \param lmaxi ...
268!> \param lmaxj ...
269!> \param irow ...
270!> \param iatom ...
271! **************************************************************************************************
272 SUBROUTINE compute_block_sk(block, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
273 llm, lmaxi, lmaxj, irow, iatom)
274 REAL(kind=dp), DIMENSION(:, :), POINTER :: block, smatij, smatji
275 REAL(kind=dp), DIMENSION(3) :: rij
276 INTEGER :: ngrd, ngrdcut
277 REAL(kind=dp) :: dgrd
278 INTEGER :: llm, lmaxi, lmaxj, irow, iatom
279
280 REAL(kind=dp) :: dr
281 REAL(kind=dp), DIMENSION(20) :: skabij, skabji
282
283 dr = sqrt(sum(rij(:)**2))
284 CALL getskz(smatij, skabij, dr, ngrd, ngrdcut, dgrd, llm)
285 CALL getskz(smatji, skabji, dr, ngrd, ngrdcut, dgrd, llm)
286 IF (irow == iatom) THEN
287 CALL turnsk(block, skabji, skabij, rij, dr, lmaxi, lmaxj)
288 ELSE
289 CALL turnsk(block, skabij, skabji, -rij, dr, lmaxj, lmaxi)
290 END IF
291
292 END SUBROUTINE compute_block_sk
293
294! **************************************************************************************************
295!> \brief Gets matrix elements on z axis, as they are stored in the tables
296!> \param slakotab ...
297!> \param skpar ...
298!> \param dx ...
299!> \param ngrd ...
300!> \param ngrdcut ...
301!> \param dgrd ...
302!> \param llm ...
303!> \author 07. Feb. 2004, TH
304! **************************************************************************************************
305 SUBROUTINE getskz(slakotab, skpar, dx, ngrd, ngrdcut, dgrd, llm)
306 REAL(dp), INTENT(in) :: slakotab(:, :), dx
307 INTEGER, INTENT(in) :: ngrd, ngrdcut
308 REAL(dp), INTENT(in) :: dgrd
309 INTEGER, INTENT(in) :: llm
310 REAL(dp), INTENT(out) :: skpar(llm)
311
312 INTEGER :: clgp
313
314 skpar = 0._dp
315 !
316 ! Determine closest grid point
317 !
318 clgp = nint(dx/dgrd)
319 !
320 ! Screen elements which are too far away
321 !
322 IF (clgp > ngrdcut) RETURN
323 !
324 ! The grid point is either contained in the table --> matrix element
325 ! can be interpolated, or it is outside the table --> matrix element
326 ! needs to be extrapolated.
327 !
328 IF (clgp > ngrd) THEN
329 !
330 ! Extrapolate external matrix elements if table does not finish with zero
331 !
332 CALL extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
333 ELSE
334 !
335 ! Interpolate tabulated matrix elements
336 !
337 CALL interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp)
338 END IF
339 END SUBROUTINE getskz
340
341! **************************************************************************************************
342!> \brief ...
343!> \param slakotab ...
344!> \param skpar ...
345!> \param dx ...
346!> \param ngrd ...
347!> \param dgrd ...
348!> \param llm ...
349!> \param clgp ...
350! **************************************************************************************************
351 SUBROUTINE interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp)
352 REAL(dp), INTENT(in) :: slakotab(:, :), dx
353 INTEGER, INTENT(in) :: ngrd
354 REAL(dp), INTENT(in) :: dgrd
355 INTEGER, INTENT(in) :: llm
356 REAL(dp), INTENT(out) :: skpar(llm)
357 INTEGER, INTENT(in) :: clgp
358
359 INTEGER :: fgpm, k, l, lgpm
360 REAL(dp) :: error, xa(max_inter), ya(max_inter)
361
362 lgpm = min(clgp + int(max_inter/2.0), ngrd)
363 fgpm = lgpm - max_inter + 1
364 DO k = 0, max_inter - 1
365 xa(k + 1) = (fgpm + k)*dgrd
366 END DO
367 !
368 ! Interpolate matrix elements for all orbitals
369 !
370 DO l = 1, llm
371 !
372 ! Read SK parameters from table
373 !
374 ya(1:max_inter) = slakotab(fgpm:lgpm, l)
375 CALL polint(xa, ya, max_inter, dx, skpar(l), error)
376 END DO
377 END SUBROUTINE interpol
378
379! **************************************************************************************************
380!> \brief ...
381!> \param slakotab ...
382!> \param skpar ...
383!> \param dx ...
384!> \param ngrd ...
385!> \param dgrd ...
386!> \param llm ...
387! **************************************************************************************************
388 SUBROUTINE extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
389 REAL(dp), INTENT(in) :: slakotab(:, :), dx
390 INTEGER, INTENT(in) :: ngrd
391 REAL(dp), INTENT(in) :: dgrd
392 INTEGER, INTENT(in) :: llm
393 REAL(dp), INTENT(out) :: skpar(llm)
394
395 INTEGER :: fgp, k, l, lgp, ntable, nzero
396 REAL(dp) :: error, xa(max_extra), ya(max_extra)
397
398 nzero = max_extra/3
399 ntable = max_extra - nzero
400 !
401 ! Get the three last distances from the table
402 !
403 DO k = 1, ntable
404 xa(k) = (ngrd - (max_extra - 3) + k)*dgrd
405 END DO
406 DO k = 1, nzero
407 xa(ntable + k) = (ngrd + k - 1)*dgrd + slako_d0
408 ya(ntable + k) = 0.0
409 END DO
410 !
411 ! Extrapolate matrix elements for all orbitals
412 !
413 DO l = 1, llm
414 !
415 ! Read SK parameters from table
416 !
417 fgp = ngrd + 1 - (max_extra - 3)
418 lgp = ngrd
419 ya(1:max_extra - 3) = slakotab(fgp:lgp, l)
420 CALL polint(xa, ya, max_extra, dx, skpar(l), error)
421 END DO
422 END SUBROUTINE extrapol
423
424! **************************************************************************************************
425!> \brief Turn matrix element from z-axis to orientation of dxv
426!> \param mat ...
427!> \param skab1 ...
428!> \param skab2 ...
429!> \param dxv ...
430!> \param dx ...
431!> \param lmaxa ...
432!> \param lmaxb ...
433!> \date 13. Jan 2004
434!> \par Notes
435!> These routines are taken from an old TB code (unknown to TH).
436!> They are highly optimised and taken because they are time critical.
437!> They are explicit, so not recursive, and work up to d functions.
438!>
439!> Set variables necessary for rotation of matrix elements
440!>
441!> r_i^2/r, replicated in rr2(4:6) for index convenience later
442!> r_i/r, direction vector, rr(4:6) are replicated from 1:3
443!> lmax of A and B
444!> \author TH
445!> \version 1.0
446! **************************************************************************************************
447 SUBROUTINE turnsk(mat, skab1, skab2, dxv, dx, lmaxa, lmaxb)
448 REAL(dp), INTENT(inout) :: mat(:, :)
449 REAL(dp), INTENT(in) :: skab1(:), skab2(:), dxv(3), dx
450 INTEGER, INTENT(in) :: lmaxa, lmaxb
451
452 INTEGER :: lmaxab, minlmaxab
453 REAL(dp) :: rinv, rr(6), rr2(6)
454
455 lmaxab = max(lmaxa, lmaxb)
456 ! Determine l quantum limits.
457 IF (lmaxab > 2) cpabort('lmax=2')
458 minlmaxab = min(lmaxa, lmaxb)
459 !
460 ! s-s interaction
461 !
462 CALL skss(skab1, mat)
463 !
464 IF (lmaxab <= 0) RETURN
465 !
466 rr2(1:3) = dxv(1:3)**2
467 rr(1:3) = dxv(1:3)
468 rinv = 1.0_dp/dx
469 !
470 rr(1:3) = rr(1:3)*rinv
471 rr(4:6) = rr(1:3)
472 rr2(1:3) = rr2(1:3)*rinv**2
473 rr2(4:6) = rr2(1:3)
474 !
475 ! s-p, p-s and p-p interaction
476 !
477 IF (minlmaxab >= 1) THEN
478 CALL skpp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
479 CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .true.)
480 CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .false.)
481 ELSE
482 IF (lmaxb >= 1) THEN
483 CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .true.)
484 ELSE
485 CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .false.)
486 END IF
487 END IF
488 !
489 ! If there is only s-p interaction we have finished
490 !
491 IF (lmaxab <= 1) RETURN
492 !
493 ! at least one atom has d functions
494 !
495 IF (minlmaxab == 2) THEN
496 !
497 ! in case both atoms have d functions
498 !
499 CALL skdd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
500 CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .true.)
501 CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .false.)
502 CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .true.)
503 CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .false.)
504 ELSE
505 !
506 ! One atom has d functions, the other has s or s and p functions
507 !
508 IF (lmaxa == 0) THEN
509 !
510 ! atom b has d, the atom a only s functions
511 !
512 CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .true.)
513 ELSE IF (lmaxa == 1) THEN
514 !
515 ! atom b has d, the atom a s and p functions
516 !
517 CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .true.)
518 CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .true.)
519 ELSE
520 !
521 ! atom a has d functions
522 !
523 IF (lmaxb == 0) THEN
524 !
525 ! atom a has d, atom b has only s functions
526 !
527 CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .false.)
528 ELSE
529 !
530 ! atom a has d, atom b has s and p functions
531 !
532 CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .false.)
533 CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .false.)
534 END IF
535 END IF
536 END IF
537 !
538 CONTAINS
539 !
540 ! The subroutines to turn the matrix elements are taken as internal subroutines
541 ! as it is beneficial to inline them.
542 !
543 ! They are both turning the matrix elements and placing them appropriately
544 ! into the matrix block
545 !
546! **************************************************************************************************
547!> \brief s-s interaction (no rotation necessary)
548!> \param skpar ...
549!> \param mat ...
550!> \version 1.0
551! **************************************************************************************************
552 SUBROUTINE skss(skpar, mat)
553 REAL(dp), INTENT(in) :: skpar(:)
554 REAL(dp), INTENT(inout) :: mat(:, :)
555
556 mat(1, 1) = mat(1, 1) + skpar(1)
557 !
558 END SUBROUTINE skss
559
560! **************************************************************************************************
561!> \brief s-p interaction (simple rotation)
562!> \param skpar ...
563!> \param mat ...
564!> \param ind ...
565!> \param transposed ...
566!> \version 1.0
567! **************************************************************************************************
568 SUBROUTINE sksp(skpar, mat, ind, transposed)
569 REAL(dp), INTENT(in) :: skpar(:)
570 REAL(dp), INTENT(inout) :: mat(:, :)
571 INTEGER, INTENT(in) :: ind(0:, 0:, 0:)
572 LOGICAL, INTENT(in) :: transposed
573
574 INTEGER :: l
575 REAL(dp) :: skp
576
577 skp = skpar(ind(1, 0, 0))
578 IF (transposed) THEN
579 DO l = 1, 3
580 mat(1, l + 1) = mat(1, l + 1) + rr(l)*skp
581 END DO
582 ELSE
583 DO l = 1, 3
584 mat(l + 1, 1) = mat(l + 1, 1) - rr(l)*skp
585 END DO
586 END IF
587 !
588 END SUBROUTINE sksp
589
590! **************************************************************************************************
591!> \brief ...
592!> \param skpar ...
593!> \param mat ...
594!> \param ind ...
595! **************************************************************************************************
596 SUBROUTINE skpp(skpar, mat, ind)
597 REAL(dp), INTENT(in) :: skpar(:)
598 REAL(dp), INTENT(inout) :: mat(:, :)
599 INTEGER, INTENT(in) :: ind(0:, 0:, 0:)
600
601 INTEGER :: ii, ir, is, k, l
602 REAL(dp) :: epp(6), matel(6), skppp, skpps
603
604 epp(1:3) = rr2(1:3)
605 DO l = 1, 3
606 epp(l + 3) = rr(l)*rr(l + 1)
607 END DO
608 skppp = skpar(ind(1, 1, 1))
609 skpps = skpar(ind(1, 1, 0))
610 !
611 DO l = 1, 3
612 matel(l) = epp(l)*skpps + (1._dp - epp(l))*skppp
613 END DO
614 DO l = 4, 6
615 matel(l) = epp(l)*(skpps - skppp)
616 END DO
617 !
618 DO ir = 1, 3
619 DO is = 1, ir - 1
620 ii = ir - is
621 k = 3*ii - (ii*(ii - 1))/2 + is
622 mat(is + 1, ir + 1) = mat(is + 1, ir + 1) + matel(k)
623 mat(ir + 1, is + 1) = mat(ir + 1, is + 1) + matel(k)
624 END DO
625 mat(ir + 1, ir + 1) = mat(ir + 1, ir + 1) + matel(ir)
626 END DO
627 END SUBROUTINE skpp
628
629! **************************************************************************************************
630!> \brief ...
631!> \param skpar ...
632!> \param mat ...
633!> \param ind ...
634!> \param transposed ...
635! **************************************************************************************************
636 SUBROUTINE sksd(skpar, mat, ind, transposed)
637 REAL(dp), INTENT(in) :: skpar(:)
638 REAL(dp), INTENT(inout) :: mat(:, :)
639 INTEGER, INTENT(in) :: ind(0:, 0:, 0:)
640 LOGICAL, INTENT(in) :: transposed
641
642 INTEGER :: l
643 REAL(dp) :: d4, d5, es(5), r3, sksds
644
645 sksds = skpar(ind(2, 0, 0))
646 r3 = sqrt(3._dp)
647 d4 = rr2(3) - 0.5_dp*(rr2(1) + rr2(2))
648 d5 = rr2(1) - rr2(2)
649 !
650 DO l = 1, 3
651 es(l) = r3*rr(l)*rr(l + 1)
652 END DO
653 es(4) = 0.5_dp*r3*d5
654 es(5) = d4
655 !
656 IF (transposed) THEN
657 DO l = 1, 5
658 mat(1, l + 4) = mat(1, l + 4) + es(l)*sksds
659 END DO
660 ELSE
661 DO l = 1, 5
662 mat(l + 4, 1) = mat(l + 4, 1) + es(l)*sksds
663 END DO
664 END IF
665 END SUBROUTINE sksd
666
667! **************************************************************************************************
668!> \brief ...
669!> \param skpar ...
670!> \param mat ...
671!> \param ind ...
672!> \param transposed ...
673! **************************************************************************************************
674 SUBROUTINE skpd(skpar, mat, ind, transposed)
675 REAL(dp), INTENT(in) :: skpar(:)
676 REAL(dp), INTENT(inout) :: mat(:, :)
677 INTEGER, INTENT(in) :: ind(0:, 0:, 0:)
678 LOGICAL, INTENT(in) :: transposed
679
680 INTEGER :: ir, is, k, l, m
681 REAL(dp) :: d3, d4, d5, d6, dm(15), epd(13, 2), r3, &
682 sktmp
683
684 r3 = sqrt(3.0_dp)
685 d3 = rr2(1) + rr2(2)
686 d4 = rr2(3) - 0.5_dp*d3
687 d5 = rr2(1) - rr2(2)
688 d6 = rr(1)*rr(2)*rr(3)
689 DO l = 1, 3
690 epd(l, 1) = r3*rr2(l)*rr(l + 1)
691 epd(l, 2) = rr(l + 1)*(1.0_dp - 2._dp*rr2(l))
692 epd(l + 4, 1) = r3*rr2(l)*rr(l + 2)
693 epd(l + 4, 2) = rr(l + 2)*(1.0_dp - 2*rr2(l))
694 epd(l + 7, 1) = 0.5_dp*r3*rr(l)*d5
695 epd(l + 10, 1) = rr(l)*d4
696 END DO
697 !
698 epd(4, 1) = r3*d6
699 epd(4, 2) = -2._dp*d6
700 epd(8, 2) = rr(1)*(1.0_dp - d5)
701 epd(9, 2) = -rr(2)*(1.0_dp + d5)
702 epd(10, 2) = -rr(3)*d5
703 epd(11, 2) = -r3*rr(1)*rr2(3)
704 epd(12, 2) = -r3*rr(2)*rr2(3)
705 epd(13, 2) = r3*rr(3)*d3
706 !
707 dm(1:15) = 0.0_dp
708 !
709 DO m = 1, 2
710 sktmp = skpar(ind(2, 1, m - 1))
711 dm(1) = dm(1) + epd(1, m)*sktmp
712 dm(2) = dm(2) + epd(6, m)*sktmp
713 dm(3) = dm(3) + epd(4, m)*sktmp
714 dm(5) = dm(5) + epd(2, m)*sktmp
715 dm(6) = dm(6) + epd(7, m)*sktmp
716 dm(7) = dm(7) + epd(5, m)*sktmp
717 dm(9) = dm(9) + epd(3, m)*sktmp
718 DO l = 8, 13
719 dm(l + 2) = dm(l + 2) + epd(l, m)*sktmp
720 END DO
721 END DO
722 !
723 dm(4) = dm(3)
724 dm(8) = dm(3)
725 !
726 IF (transposed) THEN
727 DO ir = 1, 5
728 DO is = 1, 3
729 k = 3*(ir - 1) + is
730 mat(is + 1, ir + 4) = mat(is + 1, ir + 4) + dm(k)
731 END DO
732 END DO
733 ELSE
734 DO ir = 1, 5
735 DO is = 1, 3
736 k = 3*(ir - 1) + is
737 mat(ir + 4, is + 1) = mat(ir + 4, is + 1) - dm(k)
738 END DO
739 END DO
740 END IF
741 !
742 END SUBROUTINE skpd
743
744! **************************************************************************************************
745!> \brief ...
746!> \param skpar ...
747!> \param mat ...
748!> \param ind ...
749! **************************************************************************************************
750 SUBROUTINE skdd(skpar, mat, ind)
751 REAL(dp), INTENT(in) :: skpar(:)
752 REAL(dp), INTENT(inout) :: mat(:, :)
753 INTEGER, INTENT(in) :: ind(0:, 0:, 0:)
754
755 INTEGER :: ii, ir, is, k, l, m
756 REAL(dp) :: d3, d4, d5, dd(3), dm(15), e(15, 3), r3
757
758 r3 = sqrt(3._dp)
759 d3 = rr2(1) + rr2(2)
760 d4 = rr2(3) - 0.5_dp*d3
761 d5 = rr2(1) - rr2(2)
762 DO l = 1, 3
763 e(l, 1) = rr2(l)*rr2(l + 1)
764 e(l, 2) = rr2(l) + rr2(l + 1) - 4._dp*e(l, 1)
765 e(l, 3) = rr2(l + 2) + e(l, 1)
766 e(l, 1) = 3._dp*e(l, 1)
767 END DO
768 e(4, 1) = d5**2
769 e(4, 2) = d3 - e(4, 1)
770 e(4, 3) = rr2(3) + 0.25_dp*e(4, 1)
771 e(4, 1) = 0.75_dp*e(4, 1)
772 e(5, 1) = d4**2
773 e(5, 2) = 3._dp*rr2(3)*d3
774 e(5, 3) = 0.75_dp*d3**2
775 dd(1) = rr(1)*rr(3)
776 dd(2) = rr(2)*rr(1)
777 dd(3) = rr(3)*rr(2)
778 DO l = 1, 2
779 e(l + 5, 1) = 3._dp*rr2(l + 1)*dd(l)
780 e(l + 5, 2) = dd(l)*(1._dp - 4._dp*rr2(l + 1))
781 e(l + 5, 3) = dd(l)*(rr2(l + 1) - 1._dp)
782 END DO
783 e(8, 1) = dd(1)*d5*1.5_dp
784 e(8, 2) = dd(1)*(1.0_dp - 2.0_dp*d5)
785 e(8, 3) = dd(1)*(0.5_dp*d5 - 1.0_dp)
786 e(9, 1) = d5*0.5_dp*d4*r3
787 e(9, 2) = -d5*rr2(3)*r3
788 e(9, 3) = d5*0.25_dp*(1.0_dp + rr2(3))*r3
789 e(10, 1) = rr2(1)*dd(3)*3.0_dp
790 e(10, 2) = (0.25_dp - rr2(1))*dd(3)*4.0_dp
791 e(10, 3) = dd(3)*(rr2(1) - 1.0_dp)
792 e(11, 1) = 1.5_dp*dd(3)*d5
793 e(11, 2) = -dd(3)*(1.0_dp + 2.0_dp*d5)
794 e(11, 3) = dd(3)*(1.0_dp + 0.5_dp*d5)
795 e(13, 3) = 0.5_dp*d5*dd(2)
796 e(13, 2) = -2.0_dp*dd(2)*d5
797 e(13, 1) = e(13, 3)*3.0_dp
798 e(12, 1) = d4*dd(1)*r3
799 e(14, 1) = d4*dd(3)*r3
800 e(15, 1) = d4*dd(2)*r3
801 e(15, 2) = -2.0_dp*r3*dd(2)*rr2(3)
802 e(15, 3) = 0.5_dp*r3*(1.0_dp + rr2(3))*dd(2)
803 e(14, 2) = r3*dd(3)*(d3 - rr2(3))
804 e(14, 3) = -r3*0.5_dp*dd(3)*d3
805 e(12, 2) = r3*dd(1)*(d3 - rr2(3))
806 e(12, 3) = -r3*0.5_dp*dd(1)*d3
807 !
808 dm(1:15) = 0._dp
809 DO l = 1, 15
810 DO m = 1, 3
811 dm(l) = dm(l) + e(l, m)*skpar(ind(2, 2, m - 1))
812 END DO
813 END DO
814 !
815 DO ir = 1, 5
816 DO is = 1, ir - 1
817 ii = ir - is
818 k = 5*ii - (ii*(ii - 1))/2 + is
819 mat(ir + 4, is + 4) = mat(ir + 4, is + 4) + dm(k)
820 mat(is + 4, ir + 4) = mat(is + 4, ir + 4) + dm(k)
821 END DO
822 mat(ir + 4, ir + 4) = mat(ir + 4, ir + 4) + dm(ir)
823 END DO
824 END SUBROUTINE skdd
825 !
826 END SUBROUTINE turnsk
827
828! **************************************************************************************************
829!> \brief ...
830!> \param xa ...
831!> \param ya ...
832!> \param n ...
833!> \param x ...
834!> \param y ...
835!> \param dy ...
836! **************************************************************************************************
837 SUBROUTINE polint(xa, ya, n, x, y, dy)
838 INTEGER, INTENT(in) :: n
839 REAL(dp), INTENT(in) :: ya(n), xa(n), x
840 REAL(dp), INTENT(out) :: y, dy
841
842 INTEGER :: i, m, ns
843 REAL(dp) :: c(n), d(n), den, dif, dift, ho, hp, w
844
845!
846!
847
848 ns = 1
849
850 dif = abs(x - xa(1))
851 DO i = 1, n
852 dift = abs(x - xa(i))
853 IF (dift < dif) THEN
854 ns = i
855 dif = dift
856 END IF
857 c(i) = ya(i)
858 d(i) = ya(i)
859 END DO
860 !
861 y = ya(ns)
862 ns = ns - 1
863 DO m = 1, n - 1
864 DO i = 1, n - m
865 ho = xa(i) - x
866 hp = xa(i + m) - x
867 w = c(i + 1) - d(i)
868 den = ho - hp
869 cpassert(den /= 0.0_dp)
870 den = w/den
871 d(i) = hp*den
872 c(i) = ho*den
873 END DO
874 IF (2*ns < n - m) THEN
875 dy = c(ns + 1)
876 ELSE
877 dy = d(ns)
878 ns = ns - 1
879 END IF
880 y = y + dy
881 END DO
882!
883 RETURN
884 END SUBROUTINE polint
885
886! **************************************************************************************************
887!> \brief ...
888!> \param rv ...
889!> \param r ...
890!> \param erep ...
891!> \param derep ...
892!> \param n_urpoly ...
893!> \param urep ...
894!> \param spdim ...
895!> \param s_cut ...
896!> \param srep ...
897!> \param spxr ...
898!> \param scoeff ...
899!> \param surr ...
900!> \param dograd ...
901! **************************************************************************************************
902 SUBROUTINE urep_egr(rv, r, erep, derep, &
903 n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, dograd)
904
905 REAL(dp), INTENT(in) :: rv(3), r
906 REAL(dp), INTENT(inout) :: erep, derep(3)
907 INTEGER, INTENT(in) :: n_urpoly
908 REAL(dp), INTENT(in) :: urep(:)
909 INTEGER, INTENT(in) :: spdim
910 REAL(dp), INTENT(in) :: s_cut, srep(3)
911 REAL(dp), POINTER :: spxr(:, :), scoeff(:, :)
912 REAL(dp), INTENT(in) :: surr(2)
913 LOGICAL, INTENT(in) :: dograd
914
915 INTEGER :: ic, isp, jsp, nsp
916 REAL(dp) :: de_z, rz
917
918 derep = 0._dp
919 de_z = 0._dp
920 IF (n_urpoly > 0) THEN
921 !
922 ! polynomial part
923 !
924 rz = urep(1) - r
925 IF (rz <= rtiny) RETURN
926 DO ic = 2, n_urpoly
927 erep = erep + urep(ic)*rz**(ic)
928 END DO
929 IF (dograd) THEN
930 DO ic = 2, n_urpoly
931 de_z = de_z - ic*urep(ic)*rz**(ic - 1)
932 END DO
933 END IF
934 ELSE IF (spdim > 0) THEN
935 !
936 ! spline part
937 !
938 ! This part is kind of proprietary Paderborn code and I won't reverse-engineer
939 ! everything in detail. What is obvious is documented.
940 !
941 ! This part has 4 regions:
942 ! a) very long range is screened
943 ! b) short-range is extrapolated with e-functions
944 ! ca) normal range is approximated with a spline
945 ! cb) longer range is extrapolated with an higher degree spline
946 !
947 IF (r > s_cut) RETURN ! screening (condition a)
948 !
949 IF (r < spxr(1, 1)) THEN
950 ! a) short range
951 erep = erep + exp(-srep(1)*r + srep(2)) + srep(3)
952 IF (dograd) de_z = de_z - srep(1)*exp(-srep(1)*r + srep(2))
953 ELSE
954 !
955 ! condition c). First determine between which places the spline is located:
956 !
957 ispg: DO isp = 1, spdim ! condition ca)
958 IF (r < spxr(isp, 1)) cycle ispg ! distance is smaller than this spline range
959 IF (r >= spxr(isp, 2)) cycle ispg ! distance is larger than this spline range
960 ! at this point we have found the correct spline interval
961 rz = r - spxr(isp, 1)
962 IF (isp /= spdim) THEN
963 nsp = 3 ! condition ca
964 DO jsp = 0, nsp
965 erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
966 END DO
967 IF (dograd) THEN
968 DO jsp = 1, nsp
969 de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
970 END DO
971 END IF
972 ELSE
973 nsp = 5 ! condition cb
974 DO jsp = 0, nsp
975 IF (jsp <= 3) THEN
976 erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
977 ELSE
978 erep = erep + surr(jsp - 3)*rz**(jsp)
979 END IF
980 END DO
981 IF (dograd) THEN
982 DO jsp = 1, nsp
983 IF (jsp <= 3) THEN
984 de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
985 ELSE
986 de_z = de_z + jsp*surr(jsp - 3)*rz**(jsp - 1)
987 END IF
988 END DO
989 END IF
990 END IF
991 EXIT ispg
992 END DO ispg
993 END IF
994 END IF
995 !
996 IF (dograd) THEN
997 IF (r > 1.e-12_dp) derep(1:3) = (de_z/r)*rv(1:3)
998 END IF
999
1000 END SUBROUTINE urep_egr
1001
1002END MODULE qs_dftb_utils
1003
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
qs_dftb_types
Definition of the DFTB parameter types.
Definition qs_dftb_types.F:12
qs_dftb_utils
Working with the DFTB parameter types.
Definition qs_dftb_utils.F:12
qs_dftb_utils::deallocate_dftb_atom_param
subroutine, public deallocate_dftb_atom_param(dftb_parameter)
...
Definition qs_dftb_utils.F:94
qs_dftb_utils::urep_egr
subroutine, public urep_egr(rv, r, erep, derep, n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, dograd)
...
Definition qs_dftb_utils.F:904
qs_dftb_utils::set_dftb_atom_param
subroutine, public set_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition qs_dftb_utils.F:178
qs_dftb_utils::compute_block_sk
subroutine, public compute_block_sk(block, smatij, smatji, rij, ngrd, ngrdcut, dgrd, llm, lmaxi, lmaxj, irow, iatom)
...
Definition qs_dftb_utils.F:274
qs_dftb_utils::write_dftb_atom_param
subroutine, public write_dftb_atom_param(dftb_parameter, subsys_section)
...
Definition qs_dftb_utils.F:216
qs_dftb_utils::iptr
integer, dimension(0:3, 0:3, 0:3, 0:3, 0:3), public iptr
Definition qs_dftb_utils.F:39
qs_dftb_utils::get_dftb_atom_param
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition qs_dftb_utils.F:124
qs_dftb_utils::allocate_dftb_atom_param
subroutine, public allocate_dftb_atom_param(dftb_parameter)
...
Definition qs_dftb_utils.F:62
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
qs_dftb_types::qs_dftb_atom_type
Definition qs_dftb_types.F:27