Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \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 618 : SUBROUTINE allocate_dftb_atom_param(dftb_parameter)
62 :
63 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
64 :
65 618 : IF (ASSOCIATED(dftb_parameter)) THEN
66 0 : CALL deallocate_dftb_atom_param(dftb_parameter)
67 : END IF
68 :
69 9270 : ALLOCATE (dftb_parameter)
70 :
71 618 : dftb_parameter%defined = .FALSE.
72 618 : dftb_parameter%name = ""
73 618 : dftb_parameter%typ = "NONE"
74 618 : dftb_parameter%z = -1
75 618 : dftb_parameter%zeff = -1.0_dp
76 618 : dftb_parameter%natorb = 0
77 618 : dftb_parameter%lmax = -1
78 3090 : dftb_parameter%skself = 0.0_dp
79 3090 : dftb_parameter%occupation = 0.0_dp
80 3090 : dftb_parameter%eta = 0.0_dp
81 618 : dftb_parameter%energy = 0.0_dp
82 618 : dftb_parameter%xi = 0.0_dp
83 618 : dftb_parameter%di = 0.0_dp
84 618 : dftb_parameter%rcdisp = 0.0_dp
85 618 : dftb_parameter%dudq = 0.0_dp
86 :
87 618 : END SUBROUTINE allocate_dftb_atom_param
88 :
89 : ! **************************************************************************************************
90 : !> \brief ...
91 : !> \param dftb_parameter ...
92 : ! **************************************************************************************************
93 618 : SUBROUTINE deallocate_dftb_atom_param(dftb_parameter)
94 :
95 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
96 :
97 618 : CPASSERT(ASSOCIATED(dftb_parameter))
98 618 : DEALLOCATE (dftb_parameter)
99 :
100 618 : 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 7967148 : 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 7967148 : CPASSERT(ASSOCIATED(dftb_parameter))
136 :
137 7967148 : IF (PRESENT(name)) name = dftb_parameter%name
138 7967148 : IF (PRESENT(typ)) typ = dftb_parameter%typ
139 7967148 : IF (PRESENT(defined)) defined = dftb_parameter%defined
140 7967148 : IF (PRESENT(z)) z = dftb_parameter%z
141 7967148 : IF (PRESENT(zeff)) zeff = dftb_parameter%zeff
142 7967148 : IF (PRESENT(natorb)) natorb = dftb_parameter%natorb
143 7967148 : IF (PRESENT(lmax)) lmax = dftb_parameter%lmax
144 11055336 : IF (PRESENT(skself)) skself = dftb_parameter%skself
145 37773388 : IF (PRESENT(eta)) eta = dftb_parameter%eta
146 7967148 : IF (PRESENT(energy)) energy = dftb_parameter%energy
147 7967148 : IF (PRESENT(cutoff)) cutoff = dftb_parameter%cutoff
148 7972428 : IF (PRESENT(occupation)) occupation = dftb_parameter%occupation
149 7967148 : IF (PRESENT(xi)) xi = dftb_parameter%xi
150 7967148 : IF (PRESENT(di)) di = dftb_parameter%di
151 7967148 : IF (PRESENT(rcdisp)) rcdisp = dftb_parameter%rcdisp
152 7967148 : IF (PRESENT(dudq)) dudq = dftb_parameter%dudq
153 :
154 7967148 : 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 2346 : 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 2346 : CPASSERT(ASSOCIATED(dftb_parameter))
190 :
191 2346 : IF (PRESENT(name)) dftb_parameter%name = name
192 2346 : IF (PRESENT(typ)) dftb_parameter%typ = typ
193 2346 : IF (PRESENT(defined)) dftb_parameter%defined = defined
194 2346 : IF (PRESENT(z)) dftb_parameter%z = z
195 2346 : IF (PRESENT(zeff)) dftb_parameter%zeff = zeff
196 2346 : IF (PRESENT(natorb)) dftb_parameter%natorb = natorb
197 2346 : IF (PRESENT(lmax)) dftb_parameter%lmax = lmax
198 4818 : IF (PRESENT(skself)) dftb_parameter%skself = skself
199 4818 : IF (PRESENT(eta)) dftb_parameter%eta = eta
200 4818 : IF (PRESENT(occupation)) dftb_parameter%occupation = occupation
201 2346 : IF (PRESENT(energy)) dftb_parameter%energy = energy
202 2346 : IF (PRESENT(cutoff)) dftb_parameter%cutoff = cutoff
203 2346 : IF (PRESENT(xi)) dftb_parameter%xi = xi
204 2346 : IF (PRESENT(di)) dftb_parameter%di = di
205 2346 : IF (PRESENT(rcdisp)) dftb_parameter%rcdisp = rcdisp
206 2346 : IF (PRESENT(dudq)) dftb_parameter%dudq = dudq
207 :
208 2346 : END SUBROUTINE set_dftb_atom_param
209 :
210 : ! **************************************************************************************************
211 : !> \brief ...
212 : !> \param dftb_parameter ...
213 : !> \param subsys_section ...
214 : ! **************************************************************************************************
215 618 : 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 618 : NULLIFY (logger)
227 618 : logger => cp_get_default_logger()
228 618 : 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 0 : extension=".Log")
234 :
235 0 : IF (output_unit > 0) THEN
236 : CALL get_dftb_atom_param(dftb_parameter, name=name, typ=typ, defined=defined, &
237 0 : z=z, zeff=zeff, natorb=natorb, lmax=lmax)
238 :
239 : WRITE (UNIT=output_unit, FMT="(/,A,T67,A14)") &
240 0 : " DFTB parameters: ", TRIM(name)
241 0 : IF (defined) THEN
242 : WRITE (UNIT=output_unit, FMT="(T16,A,T71,F10.2)") &
243 0 : "Effective core charge:", zeff
244 : WRITE (UNIT=output_unit, FMT="(T16,A,T71,I10)") &
245 0 : "Number of orbitals:", natorb
246 : ELSE
247 : WRITE (UNIT=output_unit, FMT="(T55,A)") &
248 0 : "Parameters are not defined"
249 : END IF
250 : END IF
251 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
252 0 : "PRINT%KINDS")
253 : END IF
254 :
255 618 : 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 6186301 : 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 24745204 : dr = SQRT(SUM(rij(:)**2))
284 6186301 : CALL getskz(smatij, skabij, dr, ngrd, ngrdcut, dgrd, llm)
285 6186301 : CALL getskz(smatji, skabji, dr, ngrd, ngrdcut, dgrd, llm)
286 6186301 : IF (irow == iatom) THEN
287 3602101 : CALL turnsk(block, skabji, skabij, rij, dr, lmaxi, lmaxj)
288 : ELSE
289 10336800 : CALL turnsk(block, skabij, skabji, -rij, dr, lmaxj, lmaxi)
290 : END IF
291 :
292 6186301 : 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 12372602 : 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 35692822 : skpar = 0._dp
315 : !
316 : ! Determine closest grid point
317 : !
318 12372602 : clgp = NINT(dx/dgrd)
319 : !
320 : ! Screen elements which are too far away
321 : !
322 12372602 : 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 12372602 : IF (clgp > ngrd) THEN
329 : !
330 : ! Extrapolate external matrix elements if table does not finish with zero
331 : !
332 2992938 : CALL extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
333 : ELSE
334 : !
335 : ! Interpolate tabulated matrix elements
336 : !
337 9379664 : 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 9379664 : 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 9379664 : lgpm = MIN(clgp + INT(max_inter/2.0), ngrd)
363 9379664 : fgpm = lgpm - max_inter + 1
364 56277984 : DO k = 0, max_inter - 1
365 56277984 : xa(k + 1) = (fgpm + k)*dgrd
366 : END DO
367 : !
368 : ! Interpolate matrix elements for all orbitals
369 : !
370 27156578 : DO l = 1, llm
371 : !
372 : ! Read SK parameters from table
373 : !
374 106661484 : ya(1:max_inter) = slakotab(fgpm:lgpm, l)
375 27156578 : CALL polint(xa, ya, max_inter, dx, skpar(l), error)
376 : END DO
377 9379664 : 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 2992938 : 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 2992938 : nzero = max_extra/3
399 2992938 : ntable = max_extra - nzero
400 : !
401 : ! Get the three last distances from the table
402 : !
403 20950566 : DO k = 1, ntable
404 20950566 : xa(k) = (ngrd - (max_extra - 3) + k)*dgrd
405 : END DO
406 11971752 : DO k = 1, nzero
407 8978814 : xa(ntable + k) = (ngrd + k - 1)*dgrd + slako_d0
408 11971752 : ya(ntable + k) = 0.0
409 : END DO
410 : !
411 : ! Extrapolate matrix elements for all orbitals
412 : !
413 8536244 : DO l = 1, llm
414 : !
415 : ! Read SK parameters from table
416 : !
417 5543306 : fgp = ngrd + 1 - (max_extra - 3)
418 5543306 : lgp = ngrd
419 38803142 : ya(1:max_extra - 3) = slakotab(fgp:lgp, l)
420 8536244 : CALL polint(xa, ya, max_extra, dx, skpar(l), error)
421 : END DO
422 2992938 : 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 6186301 : 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 6186301 : lmaxab = MAX(lmaxa, lmaxb)
456 : ! Determine l quantum limits.
457 6186301 : IF (lmaxab > 2) CPABORT('lmax=2')
458 6186301 : minlmaxab = MIN(lmaxa, lmaxb)
459 : !
460 : ! s-s interaction
461 : !
462 6186301 : CALL skss(skab1, mat)
463 : !
464 6186301 : IF (lmaxab <= 0) RETURN
465 : !
466 14176860 : rr2(1:3) = dxv(1:3)**2
467 14176860 : rr(1:3) = dxv(1:3)
468 3544215 : rinv = 1.0_dp/dx
469 : !
470 14176860 : rr(1:3) = rr(1:3)*rinv
471 14176860 : rr(4:6) = rr(1:3)
472 14176860 : rr2(1:3) = rr2(1:3)*rinv**2
473 14176860 : rr2(4:6) = rr2(1:3)
474 : !
475 : ! s-p, p-s and p-p interaction
476 : !
477 3544215 : IF (minlmaxab >= 1) THEN
478 832605 : CALL skpp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
479 832605 : CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
480 832605 : CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
481 : ELSE
482 2711610 : IF (lmaxb >= 1) THEN
483 1302670 : CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
484 : ELSE
485 1408940 : 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 3544215 : IF (lmaxab <= 1) RETURN
492 : !
493 : ! at least one atom has d functions
494 : !
495 44080 : IF (minlmaxab == 2) THEN
496 : !
497 : ! in case both atoms have d functions
498 : !
499 44048 : CALL skdd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
500 44048 : CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
501 44048 : CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
502 44048 : CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
503 44048 : 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 32 : IF (lmaxa == 0) THEN
509 : !
510 : ! atom b has d, the atom a only s functions
511 : !
512 0 : CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
513 32 : ELSE IF (lmaxa == 1) THEN
514 : !
515 : ! atom b has d, the atom a s and p functions
516 : !
517 0 : CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
518 0 : CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
519 : ELSE
520 : !
521 : ! atom a has d functions
522 : !
523 32 : IF (lmaxb == 0) THEN
524 : !
525 : ! atom a has d, atom b has only s functions
526 : !
527 0 : 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 32 : CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
533 32 : 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 6186301 : SUBROUTINE skss(skpar, mat)
553 : REAL(dp), INTENT(in) :: skpar(:)
554 : REAL(dp), INTENT(inout) :: mat(:, :)
555 :
556 6186301 : mat(1, 1) = mat(1, 1) + skpar(1)
557 : !
558 6186301 : 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 4376820 : 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 4376820 : skp = skpar(ind(1, 0, 0))
578 4376820 : IF (transposed) THEN
579 8541100 : DO l = 1, 3
580 8541100 : mat(1, l + 1) = mat(1, l + 1) + rr(l)*skp
581 : END DO
582 : ELSE
583 8966180 : DO l = 1, 3
584 8966180 : mat(l + 1, 1) = mat(l + 1, 1) - rr(l)*skp
585 : END DO
586 : END IF
587 : !
588 4376820 : END SUBROUTINE sksp
589 :
590 : ! **************************************************************************************************
591 : !> \brief ...
592 : !> \param skpar ...
593 : !> \param mat ...
594 : !> \param ind ...
595 : ! **************************************************************************************************
596 832605 : 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 3330420 : epp(1:3) = rr2(1:3)
605 3330420 : DO l = 1, 3
606 3330420 : epp(l + 3) = rr(l)*rr(l + 1)
607 : END DO
608 832605 : skppp = skpar(ind(1, 1, 1))
609 832605 : skpps = skpar(ind(1, 1, 0))
610 : !
611 3330420 : DO l = 1, 3
612 3330420 : matel(l) = epp(l)*skpps + (1._dp - epp(l))*skppp
613 : END DO
614 3330420 : DO l = 4, 6
615 3330420 : matel(l) = epp(l)*(skpps - skppp)
616 : END DO
617 : !
618 3330420 : DO ir = 1, 3
619 4995630 : DO is = 1, ir - 1
620 2497815 : ii = ir - is
621 2497815 : k = 3*ii - (ii*(ii - 1))/2 + is
622 2497815 : mat(is + 1, ir + 1) = mat(is + 1, ir + 1) + matel(k)
623 4995630 : mat(ir + 1, is + 1) = mat(ir + 1, is + 1) + matel(k)
624 : END DO
625 3330420 : mat(ir + 1, ir + 1) = mat(ir + 1, ir + 1) + matel(ir)
626 : END DO
627 832605 : END SUBROUTINE skpp
628 :
629 : ! **************************************************************************************************
630 : !> \brief ...
631 : !> \param skpar ...
632 : !> \param mat ...
633 : !> \param ind ...
634 : !> \param transposed ...
635 : ! **************************************************************************************************
636 88128 : 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 88128 : sksds = skpar(ind(2, 0, 0))
646 88128 : r3 = SQRT(3._dp)
647 88128 : d4 = rr2(3) - 0.5_dp*(rr2(1) + rr2(2))
648 88128 : d5 = rr2(1) - rr2(2)
649 : !
650 352512 : DO l = 1, 3
651 352512 : es(l) = r3*rr(l)*rr(l + 1)
652 : END DO
653 88128 : es(4) = 0.5_dp*r3*d5
654 88128 : es(5) = d4
655 : !
656 88128 : IF (transposed) THEN
657 264288 : DO l = 1, 5
658 264288 : mat(1, l + 4) = mat(1, l + 4) + es(l)*sksds
659 : END DO
660 : ELSE
661 264480 : DO l = 1, 5
662 264480 : mat(l + 4, 1) = mat(l + 4, 1) + es(l)*sksds
663 : END DO
664 : END IF
665 88128 : END SUBROUTINE sksd
666 :
667 : ! **************************************************************************************************
668 : !> \brief ...
669 : !> \param skpar ...
670 : !> \param mat ...
671 : !> \param ind ...
672 : !> \param transposed ...
673 : ! **************************************************************************************************
674 88128 : 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 88128 : r3 = SQRT(3.0_dp)
685 88128 : d3 = rr2(1) + rr2(2)
686 88128 : d4 = rr2(3) - 0.5_dp*d3
687 88128 : d5 = rr2(1) - rr2(2)
688 88128 : d6 = rr(1)*rr(2)*rr(3)
689 352512 : DO l = 1, 3
690 264384 : epd(l, 1) = r3*rr2(l)*rr(l + 1)
691 264384 : epd(l, 2) = rr(l + 1)*(1.0_dp - 2._dp*rr2(l))
692 264384 : epd(l + 4, 1) = r3*rr2(l)*rr(l + 2)
693 264384 : epd(l + 4, 2) = rr(l + 2)*(1.0_dp - 2*rr2(l))
694 264384 : epd(l + 7, 1) = 0.5_dp*r3*rr(l)*d5
695 352512 : epd(l + 10, 1) = rr(l)*d4
696 : END DO
697 : !
698 88128 : epd(4, 1) = r3*d6
699 88128 : epd(4, 2) = -2._dp*d6
700 88128 : epd(8, 2) = rr(1)*(1.0_dp - d5)
701 88128 : epd(9, 2) = -rr(2)*(1.0_dp + d5)
702 88128 : epd(10, 2) = -rr(3)*d5
703 88128 : epd(11, 2) = -r3*rr(1)*rr2(3)
704 88128 : epd(12, 2) = -r3*rr(2)*rr2(3)
705 88128 : epd(13, 2) = r3*rr(3)*d3
706 : !
707 88128 : dm(1:15) = 0.0_dp
708 : !
709 264384 : DO m = 1, 2
710 176256 : sktmp = skpar(ind(2, 1, m - 1))
711 176256 : dm(1) = dm(1) + epd(1, m)*sktmp
712 176256 : dm(2) = dm(2) + epd(6, m)*sktmp
713 176256 : dm(3) = dm(3) + epd(4, m)*sktmp
714 176256 : dm(5) = dm(5) + epd(2, m)*sktmp
715 176256 : dm(6) = dm(6) + epd(7, m)*sktmp
716 176256 : dm(7) = dm(7) + epd(5, m)*sktmp
717 176256 : dm(9) = dm(9) + epd(3, m)*sktmp
718 1321920 : DO l = 8, 13
719 1233792 : dm(l + 2) = dm(l + 2) + epd(l, m)*sktmp
720 : END DO
721 : END DO
722 : !
723 88128 : dm(4) = dm(3)
724 88128 : dm(8) = dm(3)
725 : !
726 88128 : IF (transposed) THEN
727 264288 : DO ir = 1, 5
728 925008 : DO is = 1, 3
729 660720 : k = 3*(ir - 1) + is
730 880960 : mat(is + 1, ir + 4) = mat(is + 1, ir + 4) + dm(k)
731 : END DO
732 : END DO
733 : ELSE
734 264480 : DO ir = 1, 5
735 925680 : DO is = 1, 3
736 661200 : k = 3*(ir - 1) + is
737 881600 : mat(ir + 4, is + 1) = mat(ir + 4, is + 1) - dm(k)
738 : END DO
739 : END DO
740 : END IF
741 : !
742 88128 : END SUBROUTINE skpd
743 :
744 : ! **************************************************************************************************
745 : !> \brief ...
746 : !> \param skpar ...
747 : !> \param mat ...
748 : !> \param ind ...
749 : ! **************************************************************************************************
750 44048 : 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 44048 : r3 = SQRT(3._dp)
759 44048 : d3 = rr2(1) + rr2(2)
760 44048 : d4 = rr2(3) - 0.5_dp*d3
761 44048 : d5 = rr2(1) - rr2(2)
762 176192 : DO l = 1, 3
763 132144 : e(l, 1) = rr2(l)*rr2(l + 1)
764 132144 : e(l, 2) = rr2(l) + rr2(l + 1) - 4._dp*e(l, 1)
765 132144 : e(l, 3) = rr2(l + 2) + e(l, 1)
766 176192 : e(l, 1) = 3._dp*e(l, 1)
767 : END DO
768 44048 : e(4, 1) = d5**2
769 44048 : e(4, 2) = d3 - e(4, 1)
770 44048 : e(4, 3) = rr2(3) + 0.25_dp*e(4, 1)
771 44048 : e(4, 1) = 0.75_dp*e(4, 1)
772 44048 : e(5, 1) = d4**2
773 44048 : e(5, 2) = 3._dp*rr2(3)*d3
774 44048 : e(5, 3) = 0.75_dp*d3**2
775 44048 : dd(1) = rr(1)*rr(3)
776 44048 : dd(2) = rr(2)*rr(1)
777 44048 : dd(3) = rr(3)*rr(2)
778 132144 : DO l = 1, 2
779 88096 : e(l + 5, 1) = 3._dp*rr2(l + 1)*dd(l)
780 88096 : e(l + 5, 2) = dd(l)*(1._dp - 4._dp*rr2(l + 1))
781 132144 : e(l + 5, 3) = dd(l)*(rr2(l + 1) - 1._dp)
782 : END DO
783 44048 : e(8, 1) = dd(1)*d5*1.5_dp
784 44048 : e(8, 2) = dd(1)*(1.0_dp - 2.0_dp*d5)
785 44048 : e(8, 3) = dd(1)*(0.5_dp*d5 - 1.0_dp)
786 44048 : e(9, 1) = d5*0.5_dp*d4*r3
787 44048 : e(9, 2) = -d5*rr2(3)*r3
788 44048 : e(9, 3) = d5*0.25_dp*(1.0_dp + rr2(3))*r3
789 44048 : e(10, 1) = rr2(1)*dd(3)*3.0_dp
790 44048 : e(10, 2) = (0.25_dp - rr2(1))*dd(3)*4.0_dp
791 44048 : e(10, 3) = dd(3)*(rr2(1) - 1.0_dp)
792 44048 : e(11, 1) = 1.5_dp*dd(3)*d5
793 44048 : e(11, 2) = -dd(3)*(1.0_dp + 2.0_dp*d5)
794 44048 : e(11, 3) = dd(3)*(1.0_dp + 0.5_dp*d5)
795 44048 : e(13, 3) = 0.5_dp*d5*dd(2)
796 44048 : e(13, 2) = -2.0_dp*dd(2)*d5
797 44048 : e(13, 1) = e(13, 3)*3.0_dp
798 44048 : e(12, 1) = d4*dd(1)*r3
799 44048 : e(14, 1) = d4*dd(3)*r3
800 44048 : e(15, 1) = d4*dd(2)*r3
801 44048 : e(15, 2) = -2.0_dp*r3*dd(2)*rr2(3)
802 44048 : e(15, 3) = 0.5_dp*r3*(1.0_dp + rr2(3))*dd(2)
803 44048 : e(14, 2) = r3*dd(3)*(d3 - rr2(3))
804 44048 : e(14, 3) = -r3*0.5_dp*dd(3)*d3
805 44048 : e(12, 2) = r3*dd(1)*(d3 - rr2(3))
806 44048 : e(12, 3) = -r3*0.5_dp*dd(1)*d3
807 : !
808 44048 : dm(1:15) = 0._dp
809 704768 : DO l = 1, 15
810 2686928 : DO m = 1, 3
811 2642880 : dm(l) = dm(l) + e(l, m)*skpar(ind(2, 2, m - 1))
812 : END DO
813 : END DO
814 : !
815 264288 : DO ir = 1, 5
816 660720 : DO is = 1, ir - 1
817 440480 : ii = ir - is
818 440480 : k = 5*ii - (ii*(ii - 1))/2 + is
819 440480 : mat(ir + 4, is + 4) = mat(ir + 4, is + 4) + dm(k)
820 660720 : mat(is + 4, ir + 4) = mat(is + 4, ir + 4) + dm(k)
821 : END DO
822 264288 : mat(ir + 4, ir + 4) = mat(ir + 4, ir + 4) + dm(ir)
823 : END DO
824 44048 : 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 23320220 : 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 46640440 : REAL(dp) :: c(n), d(n), den, dif, dift, ho, hp, w
844 :
845 : !
846 : !
847 :
848 23320220 : ns = 1
849 :
850 23320220 : dif = ABS(x - xa(1))
851 162094544 : DO i = 1, n
852 138774324 : dift = ABS(x - xa(i))
853 138774324 : IF (dift < dif) THEN
854 66822666 : ns = i
855 66822666 : dif = dift
856 : END IF
857 138774324 : c(i) = ya(i)
858 162094544 : d(i) = ya(i)
859 : END DO
860 : !
861 23320220 : y = ya(ns)
862 23320220 : ns = ns - 1
863 138774324 : DO m = 1, n - 1
864 492782260 : DO i = 1, n - m
865 377328156 : ho = xa(i) - x
866 377328156 : hp = xa(i + m) - x
867 377328156 : w = c(i + 1) - d(i)
868 377328156 : den = ho - hp
869 377328156 : CPASSERT(den /= 0.0_dp)
870 377328156 : den = w/den
871 377328156 : d(i) = hp*den
872 492782260 : c(i) = ho*den
873 : END DO
874 115454104 : IF (2*ns < n - m) THEN
875 48631438 : dy = c(ns + 1)
876 : ELSE
877 66822666 : dy = d(ns)
878 66822666 : ns = ns - 1
879 : END IF
880 138774324 : y = y + dy
881 : END DO
882 : !
883 23320220 : 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 581197 : SUBROUTINE urep_egr(rv, r, erep, derep, &
903 581197 : 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 581197 : derep = 0._dp
919 581197 : de_z = 0._dp
920 581197 : IF (n_urpoly > 0) THEN
921 : !
922 : ! polynomial part
923 : !
924 7761 : rz = urep(1) - r
925 7761 : IF (rz <= rtiny) RETURN
926 55144 : DO ic = 2, n_urpoly
927 55144 : erep = erep + urep(ic)*rz**(ic)
928 : END DO
929 7761 : IF (dograd) THEN
930 14087 : DO ic = 2, n_urpoly
931 14087 : de_z = de_z - ic*urep(ic)*rz**(ic - 1)
932 : END DO
933 : END IF
934 573436 : 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 573436 : IF (r > s_cut) RETURN ! screening (condition a)
948 : !
949 16896 : IF (r < spxr(1, 1)) THEN
950 : ! a) short range
951 0 : erep = erep + EXP(-srep(1)*r + srep(2)) + srep(3)
952 0 : 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 296965 : ispg: DO isp = 1, spdim ! condition ca)
958 296965 : IF (r < spxr(isp, 1)) CYCLE ispg ! distance is smaller than this spline range
959 296965 : 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 16896 : rz = r - spxr(isp, 1)
962 16896 : IF (isp /= spdim) THEN
963 : nsp = 3 ! condition ca
964 68160 : DO jsp = 0, nsp
965 68160 : erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
966 : END DO
967 13632 : IF (dograd) THEN
968 19716 : DO jsp = 1, nsp
969 19716 : 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 22848 : DO jsp = 0, nsp
975 22848 : IF (jsp <= 3) THEN
976 13056 : erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
977 : ELSE
978 6528 : erep = erep + surr(jsp - 3)*rz**(jsp)
979 : END IF
980 : END DO
981 3264 : IF (dograd) THEN
982 9792 : DO jsp = 1, nsp
983 9792 : IF (jsp <= 3) THEN
984 4896 : de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
985 : ELSE
986 3264 : de_z = de_z + jsp*surr(jsp - 3)*rz**(jsp - 1)
987 : END IF
988 : END DO
989 : END IF
990 : END IF
991 280069 : EXIT ispg
992 : END DO ispg
993 : END IF
994 : END IF
995 : !
996 24657 : IF (dograd) THEN
997 33916 : 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 :
1002 : END MODULE qs_dftb_utils
1003 :
|
