Line data Source code
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 480 : SUBROUTINE allocate_dftb_atom_param(dftb_parameter)
62 :
63 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
64 :
65 480 : IF (ASSOCIATED(dftb_parameter)) &
66 0 : CALL deallocate_dftb_atom_param(dftb_parameter)
67 :
68 7200 : ALLOCATE (dftb_parameter)
69 :
70 : dftb_parameter%defined = .FALSE.
71 480 : dftb_parameter%name = ""
72 480 : dftb_parameter%typ = "NONE"
73 : dftb_parameter%z = -1
74 : dftb_parameter%zeff = -1.0_dp
75 480 : dftb_parameter%natorb = 0
76 : dftb_parameter%lmax = -1
77 2400 : dftb_parameter%skself = 0.0_dp
78 2400 : dftb_parameter%occupation = 0.0_dp
79 2400 : dftb_parameter%eta = 0.0_dp
80 480 : dftb_parameter%energy = 0.0_dp
81 480 : dftb_parameter%xi = 0.0_dp
82 480 : dftb_parameter%di = 0.0_dp
83 480 : dftb_parameter%rcdisp = 0.0_dp
84 480 : dftb_parameter%dudq = 0.0_dp
85 :
86 480 : END SUBROUTINE allocate_dftb_atom_param
87 :
88 : ! **************************************************************************************************
89 : !> \brief ...
90 : !> \param dftb_parameter ...
91 : ! **************************************************************************************************
92 480 : SUBROUTINE deallocate_dftb_atom_param(dftb_parameter)
93 :
94 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
95 :
96 480 : CPASSERT(ASSOCIATED(dftb_parameter))
97 480 : DEALLOCATE (dftb_parameter)
98 :
99 480 : 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 6408015 : 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 6408015 : CPASSERT(ASSOCIATED(dftb_parameter))
135 :
136 6408015 : IF (PRESENT(name)) name = dftb_parameter%name
137 6408015 : IF (PRESENT(typ)) typ = dftb_parameter%typ
138 6408015 : IF (PRESENT(defined)) defined = dftb_parameter%defined
139 6408015 : IF (PRESENT(z)) z = dftb_parameter%z
140 6408015 : IF (PRESENT(zeff)) zeff = dftb_parameter%zeff
141 6408015 : IF (PRESENT(natorb)) natorb = dftb_parameter%natorb
142 6408015 : IF (PRESENT(lmax)) lmax = dftb_parameter%lmax
143 9314507 : IF (PRESENT(skself)) skself = dftb_parameter%skself
144 30314071 : IF (PRESENT(eta)) eta = dftb_parameter%eta
145 6408015 : IF (PRESENT(energy)) energy = dftb_parameter%energy
146 6408015 : IF (PRESENT(cutoff)) cutoff = dftb_parameter%cutoff
147 6410367 : IF (PRESENT(occupation)) occupation = dftb_parameter%occupation
148 6408015 : IF (PRESENT(xi)) xi = dftb_parameter%xi
149 6408015 : IF (PRESENT(di)) di = dftb_parameter%di
150 6408015 : IF (PRESENT(rcdisp)) rcdisp = dftb_parameter%rcdisp
151 6408015 : IF (PRESENT(dudq)) dudq = dftb_parameter%dudq
152 :
153 6408015 : 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 1798 : 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 1798 : CPASSERT(ASSOCIATED(dftb_parameter))
189 :
190 1798 : IF (PRESENT(name)) dftb_parameter%name = name
191 1798 : IF (PRESENT(typ)) dftb_parameter%typ = typ
192 1798 : IF (PRESENT(defined)) dftb_parameter%defined = defined
193 1798 : IF (PRESENT(z)) dftb_parameter%z = z
194 1798 : IF (PRESENT(zeff)) dftb_parameter%zeff = zeff
195 1798 : IF (PRESENT(natorb)) dftb_parameter%natorb = natorb
196 1798 : IF (PRESENT(lmax)) dftb_parameter%lmax = lmax
197 3718 : IF (PRESENT(skself)) dftb_parameter%skself = skself
198 3718 : IF (PRESENT(eta)) dftb_parameter%eta = eta
199 3718 : IF (PRESENT(occupation)) dftb_parameter%occupation = occupation
200 1798 : IF (PRESENT(energy)) dftb_parameter%energy = energy
201 1798 : IF (PRESENT(cutoff)) dftb_parameter%cutoff = cutoff
202 1798 : IF (PRESENT(xi)) dftb_parameter%xi = xi
203 1798 : IF (PRESENT(di)) dftb_parameter%di = di
204 1798 : IF (PRESENT(rcdisp)) dftb_parameter%rcdisp = rcdisp
205 1798 : IF (PRESENT(dudq)) dftb_parameter%dudq = dudq
206 :
207 1798 : END SUBROUTINE set_dftb_atom_param
208 :
209 : ! **************************************************************************************************
210 : !> \brief ...
211 : !> \param dftb_parameter ...
212 : !> \param subsys_section ...
213 : ! **************************************************************************************************
214 480 : 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 480 : NULLIFY (logger)
226 480 : logger => cp_get_default_logger()
227 480 : 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 0 : extension=".Log")
233 :
234 0 : IF (output_unit > 0) THEN
235 : CALL get_dftb_atom_param(dftb_parameter, name=name, typ=typ, defined=defined, &
236 0 : z=z, zeff=zeff, natorb=natorb, lmax=lmax)
237 :
238 : WRITE (UNIT=output_unit, FMT="(/,A,T67,A14)") &
239 0 : " DFTB parameters: ", TRIM(name)
240 0 : IF (defined) THEN
241 : WRITE (UNIT=output_unit, FMT="(T16,A,T71,F10.2)") &
242 0 : "Effective core charge:", zeff
243 : WRITE (UNIT=output_unit, FMT="(T16,A,T71,I10)") &
244 0 : "Number of orbitals:", natorb
245 : ELSE
246 : WRITE (UNIT=output_unit, FMT="(T55,A)") &
247 0 : "Parameters are not defined"
248 : END IF
249 : END IF
250 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
251 0 : "PRINT%KINDS")
252 : END IF
253 :
254 480 : 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 6036945 : 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 24147780 : dr = SQRT(SUM(rij(:)**2))
283 6036945 : CALL getskz(smatij, skabij, dr, ngrd, ngrdcut, dgrd, llm)
284 6036945 : CALL getskz(smatji, skabji, dr, ngrd, ngrdcut, dgrd, llm)
285 6036945 : IF (irow == iatom) THEN
286 3500465 : CALL turnsk(block, skabji, skabij, rij, dr, lmaxi, lmaxj)
287 : ELSE
288 10145920 : CALL turnsk(block, skabij, skabji, -rij, dr, lmaxj, lmaxi)
289 : END IF
290 :
291 6036945 : 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 12073890 : 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 34051350 : skpar = 0._dp
314 : !
315 : ! Determine closest grid point
316 : !
317 12073890 : clgp = NINT(dx/dgrd)
318 : !
319 : ! Screen elements which are too far away
320 : !
321 12073890 : 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 12072896 : IF (clgp > ngrd) THEN
328 : !
329 : ! Extrapolate external matrix elements if table does not finish with zero
330 : !
331 2922388 : CALL extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
332 : ELSE
333 : !
334 : ! Interpolate tabulated matrix elements
335 : !
336 9150508 : 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 9150508 : 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 9150508 : lgpm = MIN(clgp + INT(max_inter/2.0), ngrd)
362 9150508 : fgpm = lgpm - max_inter + 1
363 54903048 : DO k = 0, max_inter - 1
364 54903048 : xa(k + 1) = (fgpm + k)*dgrd
365 : END DO
366 : !
367 : ! Interpolate matrix elements for all orbitals
368 : !
369 25890306 : DO l = 1, llm
370 : !
371 : ! Read SK parameters from table
372 : !
373 100438788 : ya(1:max_inter) = slakotab(fgpm:lgpm, l)
374 25890306 : CALL polint(xa, ya, max_inter, dx, skpar(l), error)
375 : END DO
376 9150508 : 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 2922388 : 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 2922388 : nzero = max_extra/3
398 2922388 : ntable = max_extra - nzero
399 : !
400 : ! Get the three last distances from the table
401 : !
402 20456716 : DO k = 1, ntable
403 20456716 : xa(k) = (ngrd - (max_extra - 3) + k)*dgrd
404 : END DO
405 11689552 : DO k = 1, nzero
406 8767164 : xa(ntable + k) = (ngrd + k - 1)*dgrd + slako_d0
407 11689552 : ya(ntable + k) = 0.0
408 : END DO
409 : !
410 : ! Extrapolate matrix elements for all orbitals
411 : !
412 8158354 : DO l = 1, llm
413 : !
414 : ! Read SK parameters from table
415 : !
416 5235966 : fgp = ngrd + 1 - (max_extra - 3)
417 5235966 : lgp = ngrd
418 36651762 : ya(1:max_extra - 3) = slakotab(fgp:lgp, l)
419 8158354 : CALL polint(xa, ya, max_extra, dx, skpar(l), error)
420 : END DO
421 2922388 : 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 6036945 : 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 6036945 : lmaxab = MAX(lmaxa, lmaxb)
455 : ! Determine l quantum limits.
456 6036945 : IF (lmaxab .GT. 2) CPABORT('lmax=2')
457 6036945 : minlmaxab = MIN(lmaxa, lmaxb)
458 : !
459 : ! s-s interaction
460 : !
461 6036945 : CALL skss(skab1, mat)
462 : !
463 6036945 : IF (lmaxab .LE. 0) RETURN
464 : !
465 13648220 : rr2(1:3) = dxv(1:3)**2
466 13648220 : rr(1:3) = dxv(1:3)
467 3412055 : rinv = 1.0_dp/dx
468 : !
469 13648220 : rr(1:3) = rr(1:3)*rinv
470 13648220 : rr(4:6) = rr(1:3)
471 13648220 : rr2(1:3) = rr2(1:3)*rinv**2
472 13648220 : rr2(4:6) = rr2(1:3)
473 : !
474 : ! s-p, p-s and p-p interaction
475 : !
476 3412055 : IF (minlmaxab .GE. 1) THEN
477 714089 : CALL skpp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
478 714089 : CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
479 714089 : CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
480 : ELSE
481 2697966 : IF (lmaxb .GE. 1) THEN
482 1301534 : CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
483 : ELSE
484 1396432 : 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 3412055 : IF (lmaxab .LE. 1) RETURN
491 : !
492 : ! at least one atom has d functions
493 : !
494 18608 : IF (minlmaxab .EQ. 2) THEN
495 : !
496 : ! in case both atoms have d functions
497 : !
498 18576 : CALL skdd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
499 18576 : CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
500 18576 : CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
501 18576 : CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
502 18576 : 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 32 : IF (lmaxa .EQ. 0) THEN
508 : !
509 : ! atom b has d, the atom a only s functions
510 : !
511 0 : CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
512 32 : ELSE IF (lmaxa .EQ. 1) THEN
513 : !
514 : ! atom b has d, the atom a s and p functions
515 : !
516 0 : CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
517 0 : CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
518 : ELSE
519 : !
520 : ! atom a has d functions
521 : !
522 32 : IF (lmaxb .EQ. 0) THEN
523 : !
524 : ! atom a has d, atom b has only s functions
525 : !
526 0 : 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 32 : CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
532 32 : 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 6036945 : SUBROUTINE skss(skpar, mat)
552 : REAL(dp), INTENT(in) :: skpar(:)
553 : REAL(dp), INTENT(inout) :: mat(:, :)
554 :
555 6036945 : mat(1, 1) = mat(1, 1) + skpar(1)
556 : !
557 6036945 : 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 4126144 : 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 4126144 : skp = skpar(ind(1, 0, 0))
577 4126144 : IF (transposed) THEN
578 8062492 : DO l = 1, 3
579 8062492 : mat(1, l + 1) = mat(1, l + 1) + rr(l)*skp
580 : END DO
581 : ELSE
582 8442084 : DO l = 1, 3
583 8442084 : mat(l + 1, 1) = mat(l + 1, 1) - rr(l)*skp
584 : END DO
585 : END IF
586 : !
587 4126144 : END SUBROUTINE sksp
588 :
589 : ! **************************************************************************************************
590 : !> \brief ...
591 : !> \param skpar ...
592 : !> \param mat ...
593 : !> \param ind ...
594 : ! **************************************************************************************************
595 714089 : 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 2856356 : epp(1:3) = rr2(1:3)
604 2856356 : DO l = 1, 3
605 2856356 : epp(l + 3) = rr(l)*rr(l + 1)
606 : END DO
607 714089 : skppp = skpar(ind(1, 1, 1))
608 714089 : skpps = skpar(ind(1, 1, 0))
609 : !
610 2856356 : DO l = 1, 3
611 2856356 : matel(l) = epp(l)*skpps + (1._dp - epp(l))*skppp
612 : END DO
613 2856356 : DO l = 4, 6
614 2856356 : matel(l) = epp(l)*(skpps - skppp)
615 : END DO
616 : !
617 2856356 : DO ir = 1, 3
618 4284534 : DO is = 1, ir - 1
619 2142267 : ii = ir - is
620 2142267 : k = 3*ii - (ii*(ii - 1))/2 + is
621 2142267 : mat(is + 1, ir + 1) = mat(is + 1, ir + 1) + matel(k)
622 4284534 : mat(ir + 1, is + 1) = mat(ir + 1, is + 1) + matel(k)
623 : END DO
624 2856356 : mat(ir + 1, ir + 1) = mat(ir + 1, ir + 1) + matel(ir)
625 : END DO
626 714089 : END SUBROUTINE skpp
627 :
628 : ! **************************************************************************************************
629 : !> \brief ...
630 : !> \param skpar ...
631 : !> \param mat ...
632 : !> \param ind ...
633 : !> \param transposed ...
634 : ! **************************************************************************************************
635 37184 : 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 37184 : sksds = skpar(ind(2, 0, 0))
645 37184 : r3 = SQRT(3._dp)
646 37184 : d4 = rr2(3) - 0.5_dp*(rr2(1) + rr2(2))
647 37184 : d5 = rr2(1) - rr2(2)
648 : !
649 148736 : DO l = 1, 3
650 148736 : es(l) = r3*rr(l)*rr(l + 1)
651 : END DO
652 37184 : es(4) = 0.5_dp*r3*d5
653 37184 : es(5) = d4
654 : !
655 37184 : IF (transposed) THEN
656 111456 : DO l = 1, 5
657 111456 : mat(1, l + 4) = mat(1, l + 4) + es(l)*sksds
658 : END DO
659 : ELSE
660 111648 : DO l = 1, 5
661 111648 : mat(l + 4, 1) = mat(l + 4, 1) + es(l)*sksds
662 : END DO
663 : END IF
664 37184 : END SUBROUTINE sksd
665 :
666 : ! **************************************************************************************************
667 : !> \brief ...
668 : !> \param skpar ...
669 : !> \param mat ...
670 : !> \param ind ...
671 : !> \param transposed ...
672 : ! **************************************************************************************************
673 37184 : 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 37184 : r3 = SQRT(3.0_dp)
684 37184 : d3 = rr2(1) + rr2(2)
685 37184 : d4 = rr2(3) - 0.5_dp*d3
686 37184 : d5 = rr2(1) - rr2(2)
687 37184 : d6 = rr(1)*rr(2)*rr(3)
688 148736 : DO l = 1, 3
689 111552 : epd(l, 1) = r3*rr2(l)*rr(l + 1)
690 111552 : epd(l, 2) = rr(l + 1)*(1.0_dp - 2._dp*rr2(l))
691 111552 : epd(l + 4, 1) = r3*rr2(l)*rr(l + 2)
692 111552 : epd(l + 4, 2) = rr(l + 2)*(1.0_dp - 2*rr2(l))
693 111552 : epd(l + 7, 1) = 0.5_dp*r3*rr(l)*d5
694 148736 : epd(l + 10, 1) = rr(l)*d4
695 : END DO
696 : !
697 37184 : epd(4, 1) = r3*d6
698 37184 : epd(4, 2) = -2._dp*d6
699 37184 : epd(8, 2) = rr(1)*(1.0_dp - d5)
700 37184 : epd(9, 2) = -rr(2)*(1.0_dp + d5)
701 37184 : epd(10, 2) = -rr(3)*d5
702 37184 : epd(11, 2) = -r3*rr(1)*rr2(3)
703 37184 : epd(12, 2) = -r3*rr(2)*rr2(3)
704 37184 : epd(13, 2) = r3*rr(3)*d3
705 : !
706 37184 : dm(1:15) = 0.0_dp
707 : !
708 111552 : DO m = 1, 2
709 74368 : sktmp = skpar(ind(2, 1, m - 1))
710 74368 : dm(1) = dm(1) + epd(1, m)*sktmp
711 74368 : dm(2) = dm(2) + epd(6, m)*sktmp
712 74368 : dm(3) = dm(3) + epd(4, m)*sktmp
713 74368 : dm(5) = dm(5) + epd(2, m)*sktmp
714 74368 : dm(6) = dm(6) + epd(7, m)*sktmp
715 74368 : dm(7) = dm(7) + epd(5, m)*sktmp
716 74368 : dm(9) = dm(9) + epd(3, m)*sktmp
717 557760 : DO l = 8, 13
718 520576 : dm(l + 2) = dm(l + 2) + epd(l, m)*sktmp
719 : END DO
720 : END DO
721 : !
722 37184 : dm(4) = dm(3)
723 37184 : dm(8) = dm(3)
724 : !
725 37184 : IF (transposed) THEN
726 111456 : DO ir = 1, 5
727 390096 : DO is = 1, 3
728 278640 : k = 3*(ir - 1) + is
729 371520 : mat(is + 1, ir + 4) = mat(is + 1, ir + 4) + dm(k)
730 : END DO
731 : END DO
732 : ELSE
733 111648 : DO ir = 1, 5
734 390768 : DO is = 1, 3
735 279120 : k = 3*(ir - 1) + is
736 372160 : mat(ir + 4, is + 1) = mat(ir + 4, is + 1) - dm(k)
737 : END DO
738 : END DO
739 : END IF
740 : !
741 37184 : END SUBROUTINE skpd
742 :
743 : ! **************************************************************************************************
744 : !> \brief ...
745 : !> \param skpar ...
746 : !> \param mat ...
747 : !> \param ind ...
748 : ! **************************************************************************************************
749 18576 : 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 18576 : r3 = SQRT(3._dp)
758 18576 : d3 = rr2(1) + rr2(2)
759 18576 : d4 = rr2(3) - 0.5_dp*d3
760 18576 : d5 = rr2(1) - rr2(2)
761 74304 : DO l = 1, 3
762 55728 : e(l, 1) = rr2(l)*rr2(l + 1)
763 55728 : e(l, 2) = rr2(l) + rr2(l + 1) - 4._dp*e(l, 1)
764 55728 : e(l, 3) = rr2(l + 2) + e(l, 1)
765 74304 : e(l, 1) = 3._dp*e(l, 1)
766 : END DO
767 18576 : e(4, 1) = d5**2
768 18576 : e(4, 2) = d3 - e(4, 1)
769 18576 : e(4, 3) = rr2(3) + 0.25_dp*e(4, 1)
770 18576 : e(4, 1) = 0.75_dp*e(4, 1)
771 18576 : e(5, 1) = d4**2
772 18576 : e(5, 2) = 3._dp*rr2(3)*d3
773 18576 : e(5, 3) = 0.75_dp*d3**2
774 18576 : dd(1) = rr(1)*rr(3)
775 18576 : dd(2) = rr(2)*rr(1)
776 18576 : dd(3) = rr(3)*rr(2)
777 55728 : DO l = 1, 2
778 37152 : e(l + 5, 1) = 3._dp*rr2(l + 1)*dd(l)
779 37152 : e(l + 5, 2) = dd(l)*(1._dp - 4._dp*rr2(l + 1))
780 55728 : e(l + 5, 3) = dd(l)*(rr2(l + 1) - 1._dp)
781 : END DO
782 18576 : e(8, 1) = dd(1)*d5*1.5_dp
783 18576 : e(8, 2) = dd(1)*(1.0_dp - 2.0_dp*d5)
784 18576 : e(8, 3) = dd(1)*(0.5_dp*d5 - 1.0_dp)
785 18576 : e(9, 1) = d5*0.5_dp*d4*r3
786 18576 : e(9, 2) = -d5*rr2(3)*r3
787 18576 : e(9, 3) = d5*0.25_dp*(1.0_dp + rr2(3))*r3
788 18576 : e(10, 1) = rr2(1)*dd(3)*3.0_dp
789 18576 : e(10, 2) = (0.25_dp - rr2(1))*dd(3)*4.0_dp
790 18576 : e(10, 3) = dd(3)*(rr2(1) - 1.0_dp)
791 18576 : e(11, 1) = 1.5_dp*dd(3)*d5
792 18576 : e(11, 2) = -dd(3)*(1.0_dp + 2.0_dp*d5)
793 18576 : e(11, 3) = dd(3)*(1.0_dp + 0.5_dp*d5)
794 18576 : e(13, 3) = 0.5_dp*d5*dd(2)
795 18576 : e(13, 2) = -2.0_dp*dd(2)*d5
796 18576 : e(13, 1) = e(13, 3)*3.0_dp
797 18576 : e(12, 1) = d4*dd(1)*r3
798 18576 : e(14, 1) = d4*dd(3)*r3
799 18576 : e(15, 1) = d4*dd(2)*r3
800 18576 : e(15, 2) = -2.0_dp*r3*dd(2)*rr2(3)
801 18576 : e(15, 3) = 0.5_dp*r3*(1.0_dp + rr2(3))*dd(2)
802 18576 : e(14, 2) = r3*dd(3)*(d3 - rr2(3))
803 18576 : e(14, 3) = -r3*0.5_dp*dd(3)*d3
804 18576 : e(12, 2) = r3*dd(1)*(d3 - rr2(3))
805 18576 : e(12, 3) = -r3*0.5_dp*dd(1)*d3
806 : !
807 18576 : dm(1:15) = 0._dp
808 297216 : DO l = 1, 15
809 1133136 : DO m = 1, 3
810 1114560 : dm(l) = dm(l) + e(l, m)*skpar(ind(2, 2, m - 1))
811 : END DO
812 : END DO
813 : !
814 111456 : DO ir = 1, 5
815 278640 : DO is = 1, ir - 1
816 185760 : ii = ir - is
817 185760 : k = 5*ii - (ii*(ii - 1))/2 + is
818 185760 : mat(ir + 4, is + 4) = mat(ir + 4, is + 4) + dm(k)
819 278640 : mat(is + 4, ir + 4) = mat(is + 4, ir + 4) + dm(k)
820 : END DO
821 111456 : mat(ir + 4, ir + 4) = mat(ir + 4, ir + 4) + dm(ir)
822 : END DO
823 18576 : 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 21975764 : 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 43951528 : REAL(dp) :: c(n), d(n), den, dif, dift, ho, hp, w
843 :
844 : !
845 : !
846 :
847 21975764 : ns = 1
848 :
849 21975764 : dif = ABS(x - xa(1))
850 152798448 : DO i = 1, n
851 130822684 : dift = ABS(x - xa(i))
852 130822684 : IF (dift .LT. dif) THEN
853 62888982 : ns = i
854 62888982 : dif = dift
855 : END IF
856 130822684 : c(i) = ya(i)
857 152798448 : d(i) = ya(i)
858 : END DO
859 : !
860 21975764 : y = ya(ns)
861 21975764 : ns = ns - 1
862 130822684 : DO m = 1, n - 1
863 464739676 : DO i = 1, n - m
864 355892756 : ho = xa(i) - x
865 355892756 : hp = xa(i + m) - x
866 355892756 : w = c(i + 1) - d(i)
867 355892756 : den = ho - hp
868 355892756 : CPASSERT(den /= 0.0_dp)
869 355892756 : den = w/den
870 355892756 : d(i) = hp*den
871 464739676 : c(i) = ho*den
872 : END DO
873 108846920 : IF (2*ns .LT. n - m) THEN
874 45957938 : dy = c(ns + 1)
875 : ELSE
876 62888982 : dy = d(ns)
877 62888982 : ns = ns - 1
878 : END IF
879 130822684 : y = y + dy
880 : END DO
881 : !
882 21975764 : 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 542295 : SUBROUTINE urep_egr(rv, r, erep, derep, &
902 542295 : 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 542295 : derep = 0._dp
918 542295 : de_z = 0._dp
919 542295 : IF (n_urpoly > 0) THEN
920 : !
921 : ! polynomial part
922 : !
923 6975 : rz = urep(1) - r
924 6975 : IF (rz <= rtiny) RETURN
925 50562 : DO ic = 2, n_urpoly
926 50562 : erep = erep + urep(ic)*rz**(ic)
927 : END DO
928 6975 : IF (dograd) THEN
929 13889 : DO ic = 2, n_urpoly
930 13889 : de_z = de_z - ic*urep(ic)*rz**(ic - 1)
931 : END DO
932 : END IF
933 535320 : 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 535320 : IF (r > s_cut) RETURN ! screening (condition a)
947 : !
948 15117 : IF (r < spxr(1, 1)) THEN
949 : ! a) short range
950 0 : erep = erep + EXP(-srep(1)*r + srep(2)) + srep(3)
951 0 : 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 235687 : ispg: DO isp = 1, spdim ! condition ca)
957 235687 : IF (r < spxr(isp, 1)) CYCLE ispg ! distance is smaller than this spline range
958 235687 : 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 15117 : rz = r - spxr(isp, 1)
961 15117 : IF (isp /= spdim) THEN
962 : nsp = 3 ! condition ca
963 59265 : DO jsp = 0, nsp
964 59265 : erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
965 : END DO
966 11853 : IF (dograd) THEN
967 19344 : DO jsp = 1, nsp
968 19344 : 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 22848 : DO jsp = 0, nsp
974 22848 : IF (jsp <= 3) THEN
975 13056 : erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
976 : ELSE
977 6528 : erep = erep + surr(jsp - 3)*rz**(jsp)
978 : END IF
979 : END DO
980 3264 : IF (dograd) THEN
981 9792 : DO jsp = 1, nsp
982 9792 : IF (jsp <= 3) THEN
983 4896 : de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
984 : ELSE
985 3264 : de_z = de_z + jsp*surr(jsp - 3)*rz**(jsp - 1)
986 : END IF
987 : END DO
988 : END IF
989 : END IF
990 220570 : EXIT ispg
991 : END DO ispg
992 : END IF
993 : END IF
994 : !
995 22092 : IF (dograd) THEN
996 33408 : 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 :
|
