(git:e7e05ae)
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 ! **************************************************************************************************
12 MODULE 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 
55 CONTAINS
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 
1001 END 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