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