(git:97501a3)
Loading...
Searching...
No Matches
atom_operators.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate the atomic operator matrices
10!> \author jgh
11!> \date 03.03.2008
12!> \version 1.0
13!>
14! **************************************************************************************************
15MODULE atom_operators
16 USE ai_onecenter, ONLY: &
17 sg_coulomb, sg_erf, sg_erfc, sg_exchange, sg_gpot, sg_kinetic, sg_kinnuc, sg_nuclear, &
18 sg_overlap, sg_proj_ol, sto_kinetic, sto_nuclear, sto_overlap
19 USE atom_types, ONLY: &
20 atom_basis_gridrep, atom_basis_type, atom_compare_grids, atom_integrals, &
21 atom_potential_type, atom_state, cgto_basis, ecp_pseudo, gth_pseudo, gto_basis, lmat, &
22 no_pseudo, num_basis, release_atom_basis, sgp_pseudo, sto_basis, upf_pseudo
23 USE atom_utils, ONLY: &
24 atom_solve, contract2, contract2add, contract4, coulomb_potential_numeric, integrate_grid, &
25 numpot_matrix, slater_density, wigner_slater_functional
26 USE dkh_main, ONLY: dkh_atom_transformation
27 USE input_constants, ONLY: &
28 barrier_conf, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
29 do_sczoramp_atom, do_zoramp_atom, poly_conf
30 USE kinds, ONLY: dp
31 USE mathconstants, ONLY: gamma1,&
32 sqrt2
33 USE mathlib, ONLY: jacobi
34 USE periodic_table, ONLY: ptable
35 USE physcon, ONLY: c_light_au
36 USE qs_grid_atom, ONLY: grid_atom_type
37#include "./base/base_uses.f90"
38
39 IMPLICIT NONE
40
41 PRIVATE
42
43 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_operators'
44
45 PUBLIC :: atom_int_setup, atom_ppint_setup, atom_int_release, atom_ppint_release
46 PUBLIC :: atom_relint_setup, atom_relint_release, atom_basis_projection_overlap
47 PUBLIC :: calculate_model_potential
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief Set up atomic integrals.
53!> \param integrals atomic integrals
54!> \param basis atomic basis set
55!> \param potential pseudo-potential
56!> \param eri_coulomb setup one-centre Coulomb Electron Repulsion Integrals
57!> \param eri_exchange setup one-centre exchange Electron Repulsion Integrals
58!> \param all_nu compute integrals for all even integer parameters [0 .. 2*l]
59!> REDUNDANT, AS THIS SUBROUTINE IS NEVER INVOKED WITH all_nu = .TRUE.
60! **************************************************************************************************
61 SUBROUTINE atom_int_setup(integrals, basis, potential, &
62 eri_coulomb, eri_exchange, all_nu)
63 TYPE(atom_integrals), INTENT(INOUT) :: integrals
64 TYPE(atom_basis_type), INTENT(INOUT) :: basis
65 TYPE(atom_potential_type), INTENT(IN), OPTIONAL :: potential
66 LOGICAL, INTENT(IN), OPTIONAL :: eri_coulomb, eri_exchange, all_nu
67
68 CHARACTER(len=*), PARAMETER :: routinen = 'atom_int_setup'
69
70 INTEGER :: handle, i, ii, info, ipiv(1000), l, l1, &
71 l2, ll, lwork, m, m1, m2, mm1, mm2, n, &
72 n1, n2, nn1, nn2, nu, nx
73 REAL(kind=dp) :: om, rc, ron, sc, x
74 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot, w, work
75 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: omat, vmat
76 REAL(kind=dp), DIMENSION(:, :), POINTER :: eri
77
78 CALL timeset(routinen, handle)
79
80 IF (integrals%status == 0) THEN
81 n = maxval(basis%nbas)
82 integrals%n = basis%nbas
83
84 IF (PRESENT(eri_coulomb)) THEN
85 integrals%eri_coulomb = eri_coulomb
86 ELSE
87 integrals%eri_coulomb = .false.
88 END IF
89 IF (PRESENT(eri_exchange)) THEN
90 integrals%eri_exchange = eri_exchange
91 ELSE
92 integrals%eri_exchange = .false.
93 END IF
94 IF (PRESENT(all_nu)) THEN
95 integrals%all_nu = all_nu
96 ELSE
97 integrals%all_nu = .false.
98 END IF
99
100 NULLIFY (integrals%ovlp, integrals%kin, integrals%core, integrals%conf)
101 DO ll = 1, SIZE(integrals%ceri)
102 NULLIFY (integrals%ceri(ll)%int, integrals%eeri(ll)%int)
103 END DO
104
105 ALLOCATE (integrals%ovlp(n, n, 0:lmat))
106 integrals%ovlp = 0._dp
107
108 ALLOCATE (integrals%kin(n, n, 0:lmat))
109 integrals%kin = 0._dp
110
111 integrals%status = 1
112
113 IF (PRESENT(potential)) THEN
114 IF (potential%confinement) THEN
115 ALLOCATE (integrals%conf(n, n, 0:lmat))
116 integrals%conf = 0._dp
117 m = basis%grid%nr
118 ALLOCATE (cpot(1:m))
119 IF (potential%conf_type == poly_conf) THEN
120 rc = potential%rcon
121 sc = potential%scon
122 cpot(1:m) = (basis%grid%rad(1:m)/rc)**sc
123 ELSEIF (potential%conf_type == barrier_conf) THEN
124 om = potential%rcon
125 ron = potential%scon
126 rc = ron + om
127 DO i = 1, m
128 IF (basis%grid%rad(i) < ron) THEN
129 cpot(i) = 0.0_dp
130 ELSEIF (basis%grid%rad(i) < rc) THEN
131 x = (basis%grid%rad(i) - ron)/om
132 x = 1._dp - x
133 cpot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
134 x = (rc - basis%grid%rad(i))**2/om/(basis%grid%rad(i) - ron)
135 cpot(i) = cpot(i)*x
136 ELSE
137 cpot(i) = 1.0_dp
138 END IF
139 END DO
140 ELSE
141 cpabort("")
142 END IF
143 CALL numpot_matrix(integrals%conf, cpot, basis, 0)
144 DEALLOCATE (cpot)
145 END IF
146 END IF
147
148 SELECT CASE (basis%basis_type)
149 CASE DEFAULT
150 cpabort("")
151 CASE (gto_basis)
152 DO l = 0, lmat
153 n = integrals%n(l)
154 CALL sg_overlap(integrals%ovlp(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
155 CALL sg_kinetic(integrals%kin(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
156 END DO
157 IF (integrals%eri_coulomb) THEN
158 ll = 0
159 DO l1 = 0, lmat
160 n1 = integrals%n(l1)
161 nn1 = (n1*(n1 + 1))/2
162 DO l2 = 0, l1
163 n2 = integrals%n(l2)
164 nn2 = (n2*(n2 + 1))/2
165 IF (integrals%all_nu) THEN
166 nx = min(2*l1, 2*l2)
167 ELSE
168 nx = 0
169 END IF
170 DO nu = 0, nx, 2
171 ll = ll + 1
172 cpassert(ll <= SIZE(integrals%ceri))
173 ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
174 integrals%ceri(ll)%int = 0._dp
175 eri => integrals%ceri(ll)%int
176 CALL sg_coulomb(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
177 END DO
178 END DO
179 END DO
180 END IF
181 IF (integrals%eri_exchange) THEN
182 ll = 0
183 DO l1 = 0, lmat
184 n1 = integrals%n(l1)
185 nn1 = (n1*(n1 + 1))/2
186 DO l2 = 0, l1
187 n2 = integrals%n(l2)
188 nn2 = (n2*(n2 + 1))/2
189 DO nu = abs(l1 - l2), l1 + l2, 2
190 ll = ll + 1
191 cpassert(ll <= SIZE(integrals%eeri))
192 ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
193 integrals%eeri(ll)%int = 0._dp
194 eri => integrals%eeri(ll)%int
195 CALL sg_exchange(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
196 END DO
197 END DO
198 END DO
199 END IF
200 CASE (cgto_basis)
201 DO l = 0, lmat
202 n = integrals%n(l)
203 m = basis%nprim(l)
204 IF (n > 0 .AND. m > 0) THEN
205 ALLOCATE (omat(m, m))
206
207 CALL sg_overlap(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
208 CALL contract2(integrals%ovlp(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
209 CALL sg_kinetic(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
210 CALL contract2(integrals%kin(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
211 DEALLOCATE (omat)
212 END IF
213 END DO
214 IF (integrals%eri_coulomb) THEN
215 ll = 0
216 DO l1 = 0, lmat
217 n1 = integrals%n(l1)
218 nn1 = (n1*(n1 + 1))/2
219 m1 = basis%nprim(l1)
220 mm1 = (m1*(m1 + 1))/2
221 DO l2 = 0, l1
222 n2 = integrals%n(l2)
223 nn2 = (n2*(n2 + 1))/2
224 m2 = basis%nprim(l2)
225 mm2 = (m2*(m2 + 1))/2
226 IF (integrals%all_nu) THEN
227 nx = min(2*l1, 2*l2)
228 ELSE
229 nx = 0
230 END IF
231 DO nu = 0, nx, 2
232 ll = ll + 1
233 cpassert(ll <= SIZE(integrals%ceri))
234 ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
235 integrals%ceri(ll)%int = 0._dp
236 ALLOCATE (omat(mm1, mm2))
237
238 eri => integrals%ceri(ll)%int
239 CALL sg_coulomb(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
240 CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))
241
242 DEALLOCATE (omat)
243 END DO
244 END DO
245 END DO
246 END IF
247 IF (integrals%eri_exchange) THEN
248 ll = 0
249 DO l1 = 0, lmat
250 n1 = integrals%n(l1)
251 nn1 = (n1*(n1 + 1))/2
252 m1 = basis%nprim(l1)
253 mm1 = (m1*(m1 + 1))/2
254 DO l2 = 0, l1
255 n2 = integrals%n(l2)
256 nn2 = (n2*(n2 + 1))/2
257 m2 = basis%nprim(l2)
258 mm2 = (m2*(m2 + 1))/2
259 DO nu = abs(l1 - l2), l1 + l2, 2
260 ll = ll + 1
261 cpassert(ll <= SIZE(integrals%eeri))
262 ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
263 integrals%eeri(ll)%int = 0._dp
264 ALLOCATE (omat(mm1, mm2))
265
266 eri => integrals%eeri(ll)%int
267 CALL sg_exchange(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
268 CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))
269
270 DEALLOCATE (omat)
271 END DO
272 END DO
273 END DO
274 END IF
275 CASE (sto_basis)
276 DO l = 0, lmat
277 n = integrals%n(l)
278 CALL sto_overlap(integrals%ovlp(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
279 basis%ns(1:n, l), basis%as(1:n, l))
280 CALL sto_kinetic(integrals%kin(1:n, 1:n, l), l, basis%ns(1:n, l), basis%as(1:n, l), &
281 basis%ns(1:n, l), basis%as(1:n, l))
282 END DO
283 cpassert(.NOT. integrals%eri_coulomb)
284 cpassert(.NOT. integrals%eri_exchange)
285 CASE (num_basis)
286 cpabort("")
287 END SELECT
288
289 ! setup transformation matrix to get an orthogonal basis, remove linear dependencies
290 NULLIFY (integrals%utrans, integrals%uptrans)
291 n = maxval(basis%nbas)
292 ALLOCATE (integrals%utrans(n, n, 0:lmat), integrals%uptrans(n, n, 0:lmat))
293 integrals%utrans = 0._dp
294 integrals%uptrans = 0._dp
295 integrals%nne = integrals%n
296 lwork = 10*n
297 ALLOCATE (omat(n, n), vmat(n, n), w(n), work(lwork))
298 DO l = 0, lmat
299 n = integrals%n(l)
300 IF (n > 0) THEN
301 omat(1:n, 1:n) = integrals%ovlp(1:n, 1:n, l)
302 CALL jacobi(omat(1:n, 1:n), w(1:n), vmat(1:n, 1:n))
303 omat(1:n, 1:n) = vmat(1:n, 1:n)
304 ii = 0
305 DO i = 1, n
306 IF (w(i) > basis%eps_eig) THEN
307 ii = ii + 1
308 integrals%utrans(1:n, ii, l) = omat(1:n, i)/sqrt(w(i))
309 END IF
310 END DO
311 integrals%nne(l) = ii
312 IF (ii > 0) THEN
313 omat(1:ii, 1:ii) = matmul(transpose(integrals%utrans(1:n, 1:ii, l)), integrals%utrans(1:n, 1:ii, l))
314 DO i = 1, ii
315 integrals%uptrans(i, i, l) = 1._dp
316 END DO
317 CALL dgesv(ii, ii, omat(1:ii, 1:ii), ii, ipiv, integrals%uptrans(1:ii, 1:ii, l), ii, info)
318 cpassert(info == 0)
319 END IF
320 END IF
321 END DO
322 DEALLOCATE (omat, vmat, w, work)
323 END IF
324
325 CALL timestop(handle)
326
327 END SUBROUTINE atom_int_setup
328
329! **************************************************************************************************
330!> \brief ...
331!> \param integrals ...
332!> \param basis ...
333!> \param potential ...
334! **************************************************************************************************
335 SUBROUTINE atom_ppint_setup(integrals, basis, potential)
336 TYPE(atom_integrals), INTENT(INOUT) :: integrals
337 TYPE(atom_basis_type), INTENT(INOUT) :: basis
338 TYPE(atom_potential_type), INTENT(IN) :: potential
339
340 CHARACTER(len=*), PARAMETER :: routinen = 'atom_ppint_setup'
341
342 INTEGER :: handle, i, ii, j, k, l, m, n
343 REAL(kind=dp) :: al, alpha, rc
344 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot, xmat
345 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: omat, spmat
346 REAL(kind=dp), DIMENSION(:), POINTER :: rad
347
348 CALL timeset(routinen, handle)
349
350 IF (integrals%ppstat == 0) THEN
351 n = maxval(basis%nbas)
352 integrals%n = basis%nbas
353
354 NULLIFY (integrals%core, integrals%hnl)
355
356 ALLOCATE (integrals%hnl(n, n, 0:lmat))
357 integrals%hnl = 0._dp
358
359 ALLOCATE (integrals%core(n, n, 0:lmat))
360 integrals%core = 0._dp
361
362 ALLOCATE (integrals%clsd(n, n, 0:lmat))
363 integrals%clsd = 0._dp
364
365 integrals%ppstat = 1
366
367 SELECT CASE (basis%basis_type)
368 CASE DEFAULT
369 cpabort("")
370 CASE (gto_basis)
371
372 SELECT CASE (potential%ppot_type)
373 CASE (no_pseudo, ecp_pseudo)
374 DO l = 0, lmat
375 n = integrals%n(l)
376 CALL sg_nuclear(integrals%core(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
377 END DO
378 CASE (gth_pseudo)
379 IF (potential%gth_pot%soc) THEN
380 cpwarn("Atom code: GTH SOC is ignored")
381 END IF
382 alpha = 1._dp/potential%gth_pot%rc/sqrt(2._dp)
383 DO l = 0, lmat
384 n = integrals%n(l)
385 ALLOCATE (omat(n, n), spmat(n, 5))
386
387 omat = 0._dp
388 CALL sg_erf(omat(1:n, 1:n), l, alpha, basis%am(1:n, l), basis%am(1:n, l))
389 integrals%core(1:n, 1:n, l) = -potential%gth_pot%zion*omat(1:n, 1:n)
390 DO i = 1, potential%gth_pot%ncl
391 omat = 0._dp
392 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%rc, l, basis%am(1:n, l), basis%am(1:n, l))
393 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
394 potential%gth_pot%cl(i)*omat(1:n, 1:n)
395 END DO
396 IF (potential%gth_pot%lpotextended) THEN
397 DO k = 1, potential%gth_pot%nexp_lpot
398 DO i = 1, potential%gth_pot%nct_lpot(k)
399 omat = 0._dp
400 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lpot(k), l, &
401 basis%am(1:n, l), basis%am(1:n, l))
402 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
403 potential%gth_pot%cval_lpot(i, k)*omat(1:n, 1:n)
404 END DO
405 END DO
406 END IF
407 IF (potential%gth_pot%lsdpot) THEN
408 DO k = 1, potential%gth_pot%nexp_lsd
409 DO i = 1, potential%gth_pot%nct_lsd(k)
410 omat = 0._dp
411 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lsd(k), l, &
412 basis%am(1:n, l), basis%am(1:n, l))
413 integrals%clsd(1:n, 1:n, l) = integrals%clsd(1:n, 1:n, l) + &
414 potential%gth_pot%cval_lsd(i, k)*omat(1:n, 1:n)
415 END DO
416 END DO
417 END IF
418
419 spmat = 0._dp
420 m = potential%gth_pot%nl(l)
421 DO i = 1, m
422 CALL sg_proj_ol(spmat(1:n, i), l, basis%am(1:n, l), i - 1, potential%gth_pot%rcnl(l))
423 END DO
424 integrals%hnl(1:n, 1:n, l) = matmul(spmat(1:n, 1:m), &
425 matmul(potential%gth_pot%hnl(1:m, 1:m, l), transpose(spmat(1:n, 1:m))))
426
427 DEALLOCATE (omat, spmat)
428 END DO
429 CASE (upf_pseudo)
430 CALL upfint_setup(integrals, basis, potential)
431 CASE (sgp_pseudo)
432 CALL sgpint_setup(integrals, basis, potential)
433 CASE DEFAULT
434 cpabort("")
435 END SELECT
436
437 CASE (cgto_basis)
438
439 SELECT CASE (potential%ppot_type)
440 CASE (no_pseudo, ecp_pseudo)
441 DO l = 0, lmat
442 n = integrals%n(l)
443 m = basis%nprim(l)
444 ALLOCATE (omat(m, m))
445
446 CALL sg_nuclear(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
447 CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
448
449 DEALLOCATE (omat)
450 END DO
451 CASE (gth_pseudo)
452 alpha = 1._dp/potential%gth_pot%rc/sqrt(2._dp)
453 DO l = 0, lmat
454 n = integrals%n(l)
455 m = basis%nprim(l)
456 IF (n > 0 .AND. m > 0) THEN
457 ALLOCATE (omat(m, m), spmat(n, 5), xmat(m))
458
459 omat = 0._dp
460 CALL sg_erf(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
461 omat(1:m, 1:m) = -potential%gth_pot%zion*omat(1:m, 1:m)
462 CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
463 DO i = 1, potential%gth_pot%ncl
464 omat = 0._dp
465 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%rc, l, basis%am(1:m, l), basis%am(1:m, l))
466 omat(1:m, 1:m) = potential%gth_pot%cl(i)*omat(1:m, 1:m)
467 CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
468 END DO
469 IF (potential%gth_pot%lpotextended) THEN
470 DO k = 1, potential%gth_pot%nexp_lpot
471 DO i = 1, potential%gth_pot%nct_lpot(k)
472 omat = 0._dp
473 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lpot(k), l, &
474 basis%am(1:m, l), basis%am(1:m, l))
475 omat(1:m, 1:m) = potential%gth_pot%cval_lpot(i, k)*omat(1:m, 1:m)
476 CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
477 END DO
478 END DO
479 END IF
480 IF (potential%gth_pot%lsdpot) THEN
481 DO k = 1, potential%gth_pot%nexp_lsd
482 DO i = 1, potential%gth_pot%nct_lsd(k)
483 omat = 0._dp
484 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lsd(k), l, &
485 basis%am(1:m, l), basis%am(1:m, l))
486 omat(1:m, 1:m) = potential%gth_pot%cval_lsd(i, k)*omat(1:m, 1:m)
487 CALL contract2add(integrals%clsd(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
488 END DO
489 END DO
490 END IF
491
492 spmat = 0._dp
493 k = potential%gth_pot%nl(l)
494 DO i = 1, k
495 CALL sg_proj_ol(xmat(1:m), l, basis%am(1:m, l), i - 1, potential%gth_pot%rcnl(l))
496 spmat(1:n, i) = matmul(transpose(basis%cm(1:m, 1:n, l)), xmat(1:m))
497 END DO
498 IF (k > 0) THEN
499 integrals%hnl(1:n, 1:n, l) = matmul(spmat(1:n, 1:k), &
500 matmul(potential%gth_pot%hnl(1:k, 1:k, l), &
501 transpose(spmat(1:n, 1:k))))
502 END IF
503 DEALLOCATE (omat, spmat, xmat)
504 END IF
505 END DO
506 CASE (upf_pseudo)
507 CALL upfint_setup(integrals, basis, potential)
508 CASE (sgp_pseudo)
509 CALL sgpint_setup(integrals, basis, potential)
510 CASE DEFAULT
511 cpabort("")
512 END SELECT
513
514 CASE (sto_basis)
515
516 SELECT CASE (potential%ppot_type)
517 CASE (no_pseudo, ecp_pseudo)
518 DO l = 0, lmat
519 n = integrals%n(l)
520 CALL sto_nuclear(integrals%core(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
521 basis%ns(1:n, l), basis%as(1:n, l))
522 END DO
523 CASE (gth_pseudo)
524 rad => basis%grid%rad
525 m = basis%grid%nr
526 ALLOCATE (cpot(1:m))
527 rc = potential%gth_pot%rc
528 alpha = 1._dp/rc/sqrt(2._dp)
529
530 ! local pseudopotential, we use erf = 1/r - erfc
531 integrals%core = 0._dp
532 DO i = 1, m
533 cpot(i) = potential%gth_pot%zion*erfc(alpha*rad(i))/rad(i)
534 END DO
535 DO i = 1, potential%gth_pot%ncl
536 ii = 2*(i - 1)
537 cpot(1:m) = cpot(1:m) + potential%gth_pot%cl(i)*(rad/rc)**ii*exp(-0.5_dp*(rad/rc)**2)
538 END DO
539 IF (potential%gth_pot%lpotextended) THEN
540 DO k = 1, potential%gth_pot%nexp_lpot
541 al = potential%gth_pot%alpha_lpot(k)
542 DO i = 1, potential%gth_pot%nct_lpot(k)
543 ii = 2*(i - 1)
544 cpot(1:m) = cpot(1:m) + potential%gth_pot%cval_lpot(i, k)*(rad/al)**ii*exp(-0.5_dp*(rad/al)**2)
545 END DO
546 END DO
547 END IF
548 CALL numpot_matrix(integrals%core, cpot, basis, 0)
549 DO l = 0, lmat
550 n = integrals%n(l)
551 ALLOCATE (omat(n, n))
552 omat = 0._dp
553 CALL sto_nuclear(omat(1:n, 1:n), basis%ns(1:n, l), basis%as(1:n, l), &
554 basis%ns(1:n, l), basis%as(1:n, l))
555 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) - potential%gth_pot%zion*omat(1:n, 1:n)
556 DEALLOCATE (omat)
557 END DO
558
559 IF (potential%gth_pot%lsdpot) THEN
560 cpot = 0._dp
561 DO k = 1, potential%gth_pot%nexp_lsd
562 al = potential%gth_pot%alpha_lsd(k)
563 DO i = 1, potential%gth_pot%nct_lsd(k)
564 ii = 2*(i - 1)
565 cpot(:) = cpot + potential%gth_pot%cval_lsd(i, k)*(rad/al)**ii*exp(-0.5_dp*(rad/al)**2)
566 END DO
567 END DO
568 CALL numpot_matrix(integrals%clsd, cpot, basis, 0)
569 END IF
570
571 DO l = 0, lmat
572 n = integrals%n(l)
573 ! non local pseudopotential
574 ALLOCATE (spmat(n, 5))
575 spmat = 0._dp
576 k = potential%gth_pot%nl(l)
577 DO i = 1, k
578 rc = potential%gth_pot%rcnl(l)
579 cpot(:) = sqrt2/sqrt(gamma1(l + 2*i - 1))*rad**(l + 2*i - 2)*exp(-0.5_dp*(rad/rc)**2)/rc**(l + 2*i - 0.5_dp)
580 DO j = 1, basis%nbas(l)
581 spmat(j, i) = integrate_grid(cpot, basis%bf(:, j, l), basis%grid)
582 END DO
583 END DO
584 integrals%hnl(1:n, 1:n, l) = matmul(spmat(1:n, 1:k), &
585 matmul(potential%gth_pot%hnl(1:k, 1:k, l), &
586 transpose(spmat(1:n, 1:k))))
587 DEALLOCATE (spmat)
588 END DO
589 DEALLOCATE (cpot)
590
591 CASE (upf_pseudo)
592 CALL upfint_setup(integrals, basis, potential)
593 CASE (sgp_pseudo)
594 CALL sgpint_setup(integrals, basis, potential)
595 CASE DEFAULT
596 cpabort("")
597 END SELECT
598
599 CASE (num_basis)
600 cpabort("")
601 END SELECT
602
603 ! add ecp_pseudo using numerical representation of basis
604 IF (potential%ppot_type == ecp_pseudo) THEN
605 ! scale 1/r potential
606 IF (any(potential%ecp_pot%nrloc(1:potential%ecp_pot%nloc) == 1)) THEN
607 alpha = 0.0_dp
608 DO k = 1, potential%ecp_pot%nloc
609 n = potential%ecp_pot%nrloc(k)
610 IF (n == 1) alpha = alpha + potential%ecp_pot%aloc(k)
611 END DO
612 integrals%core = (alpha - potential%ecp_pot%zion)*integrals%core
613 ELSE
614 integrals%core = -potential%ecp_pot%zion*integrals%core
615 END IF
616 ! local potential
617 m = basis%grid%nr
618 rad => basis%grid%rad
619 ALLOCATE (cpot(1:m))
620 cpot = 0._dp
621 DO k = 1, potential%ecp_pot%nloc
622 n = potential%ecp_pot%nrloc(k)
623 alpha = potential%ecp_pot%bloc(k)
624 IF (n == 1) THEN
625 cpot(:) = cpot + potential%ecp_pot%aloc(k)/rad*(exp(-alpha*rad**2) - 1.0_dp)
626 ELSE
627 cpot(:) = cpot + potential%ecp_pot%aloc(k)*rad**(n - 2)*exp(-alpha*rad**2)
628 END IF
629 END DO
630 CALL numpot_matrix(integrals%core, cpot, basis, 0)
631 ! non local pseudopotential
632 DO l = 0, min(potential%ecp_pot%lmax, lmat)
633 cpot = 0._dp
634 DO k = 1, potential%ecp_pot%npot(l)
635 n = potential%ecp_pot%nrpot(k, l)
636 alpha = potential%ecp_pot%bpot(k, l)
637 cpot(:) = cpot + potential%ecp_pot%apot(k, l)*rad**(n - 2)*exp(-alpha*rad**2)
638 END DO
639 DO i = 1, basis%nbas(l)
640 DO j = i, basis%nbas(l)
641 integrals%hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
642 integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
643 END DO
644 END DO
645 END DO
646 DEALLOCATE (cpot)
647 END IF
648
649 END IF
650
651 CALL timestop(handle)
652
653 END SUBROUTINE atom_ppint_setup
654
655! **************************************************************************************************
656!> \brief ...
657!> \param integrals ...
658!> \param basis ...
659!> \param potential ...
660! **************************************************************************************************
661 SUBROUTINE upfint_setup(integrals, basis, potential)
662 TYPE(atom_integrals), INTENT(INOUT) :: integrals
663 TYPE(atom_basis_type), INTENT(INOUT) :: basis
664 TYPE(atom_potential_type), INTENT(IN) :: potential
665
666 CHARACTER(len=4) :: ptype
667 INTEGER :: i, j, k1, k2, la, lb, m, n
668 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: spot
669 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: spmat
670 TYPE(atom_basis_type) :: gbasis
671
672 ! get basis representation on UPF grid
673 CALL atom_basis_gridrep(basis, gbasis, potential%upf_pot%r, potential%upf_pot%rab)
674
675 ! local pseudopotential
676 integrals%core = 0._dp
677 CALL numpot_matrix(integrals%core, potential%upf_pot%vlocal, gbasis, 0)
678
679 ptype = adjustl(trim(potential%upf_pot%pseudo_type))
680 IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
681 ! non local pseudopotential
682 n = maxval(integrals%n(:))
683 m = potential%upf_pot%number_of_proj
684 ALLOCATE (spmat(n, m))
685 spmat = 0.0_dp
686 DO i = 1, m
687 la = potential%upf_pot%lbeta(i)
688 DO j = 1, gbasis%nbas(la)
689 spmat(j, i) = integrate_grid(potential%upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
690 END DO
691 END DO
692 DO i = 1, m
693 la = potential%upf_pot%lbeta(i)
694 DO j = 1, m
695 lb = potential%upf_pot%lbeta(j)
696 IF (la == lb) THEN
697 DO k1 = 1, gbasis%nbas(la)
698 DO k2 = 1, gbasis%nbas(la)
699 integrals%hnl(k1, k2, la) = integrals%hnl(k1, k2, la) + &
700 spmat(k1, i)*potential%upf_pot%dion(i, j)*spmat(k2, j)
701 END DO
702 END DO
703 END IF
704 END DO
705 END DO
706 DEALLOCATE (spmat)
707 ELSE IF (ptype(1:2) == "SL") THEN
708 ! semi local pseudopotential
709 DO la = 0, potential%upf_pot%l_max
710 IF (la == potential%upf_pot%l_local) cycle
711 m = SIZE(potential%upf_pot%vsemi(:, la + 1))
712 ALLOCATE (spot(m))
713 spot(:) = potential%upf_pot%vsemi(:, la + 1) - potential%upf_pot%vlocal(:)
714 n = basis%nbas(la)
715 DO i = 1, n
716 DO j = i, n
717 integrals%core(i, j, la) = integrals%core(i, j, la) + &
718 integrate_grid(spot(:), &
719 gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
720 integrals%core(j, i, la) = integrals%core(i, j, la)
721 END DO
722 END DO
723 DEALLOCATE (spot)
724 END DO
725 ELSE
726 cpabort("Pseudopotential type: ["//adjustl(trim(ptype))//"] not known")
727 END IF
728
729 ! release basis representation on UPF grid
730 CALL release_atom_basis(gbasis)
731
732 END SUBROUTINE upfint_setup
733
734! **************************************************************************************************
735!> \brief ...
736!> \param integrals ...
737!> \param basis ...
738!> \param potential ...
739! **************************************************************************************************
740 SUBROUTINE sgpint_setup(integrals, basis, potential)
741 TYPE(atom_integrals), INTENT(INOUT) :: integrals
742 TYPE(atom_basis_type), INTENT(INOUT) :: basis
743 TYPE(atom_potential_type), INTENT(IN) :: potential
744
745 INTEGER :: i, ia, j, l, m, n, na
746 REAL(kind=dp) :: a, c, rc, zval
747 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot, pgauss
748 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: qmat
749 REAL(kind=dp), DIMENSION(:), POINTER :: rad
750
751 rad => basis%grid%rad
752 m = basis%grid%nr
753
754 ! local pseudopotential
755 integrals%core = 0._dp
756 ALLOCATE (cpot(m))
757 cpot = 0.0_dp
758 zval = potential%sgp_pot%zion
759 DO i = 1, m
760 rc = rad(i)/potential%sgp_pot%ac_local/sqrt(2.0_dp)
761 cpot(i) = cpot(i) - zval/rad(i)*erf(rc)
762 END DO
763 DO i = 1, potential%sgp_pot%n_local
764 cpot(:) = cpot(:) + potential%sgp_pot%c_local(i)*exp(-potential%sgp_pot%a_local(i)*rad(:)**2)
765 END DO
766 CALL numpot_matrix(integrals%core, cpot, basis, 0)
767 DEALLOCATE (cpot)
768
769 ! nonlocal pseudopotential
770 integrals%hnl = 0.0_dp
771 IF (potential%sgp_pot%has_nonlocal) THEN
772 ALLOCATE (pgauss(1:m))
773 n = potential%sgp_pot%n_nonlocal
774 !
775 DO l = 0, potential%sgp_pot%lmax
776 cpassert(l <= ubound(basis%nbas, 1))
777 IF (.NOT. potential%sgp_pot%is_nonlocal(l)) cycle
778 ! overlap (a|p)
779 na = basis%nbas(l)
780 ALLOCATE (qmat(na, n))
781 DO i = 1, n
782 pgauss(:) = 0.0_dp
783 DO j = 1, n
784 a = potential%sgp_pot%a_nonlocal(j)
785 c = potential%sgp_pot%c_nonlocal(j, i, l)
786 pgauss(:) = pgauss(:) + c*exp(-a*rad(:)**2)*rad(:)**l
787 END DO
788 DO ia = 1, na
789 qmat(ia, i) = sum(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:))
790 END DO
791 END DO
792 DO i = 1, na
793 DO j = i, na
794 DO ia = 1, n
795 integrals%hnl(i, j, l) = integrals%hnl(i, j, l) &
796 + qmat(i, ia)*qmat(j, ia)*potential%sgp_pot%h_nonlocal(ia, l)
797 END DO
798 integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
799 END DO
800 END DO
801 DEALLOCATE (qmat)
802 END DO
803 DEALLOCATE (pgauss)
804 END IF
805
806 END SUBROUTINE sgpint_setup
807
808! **************************************************************************************************
809!> \brief ...
810!> \param integrals ...
811!> \param basis ...
812!> \param reltyp ...
813!> \param zcore ...
814!> \param alpha ...
815! **************************************************************************************************
816 SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
817 TYPE(atom_integrals), INTENT(INOUT) :: integrals
818 TYPE(atom_basis_type), INTENT(INOUT) :: basis
819 INTEGER, INTENT(IN) :: reltyp
820 REAL(dp), INTENT(IN) :: zcore
821 REAL(dp), INTENT(IN), OPTIONAL :: alpha
822
823 CHARACTER(len=*), PARAMETER :: routinen = 'atom_relint_setup'
824
825 INTEGER :: dkhorder, handle, i, k1, k2, l, m, n, nl
826 REAL(dp) :: ascal
827 REAL(dp), ALLOCATABLE, DIMENSION(:) :: cpot
828 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: modpot
829 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ener, sps
830 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hmat, pvp, sp, tp, vp, wfn
831
832 CALL timeset(routinen, handle)
833
834 SELECT CASE (reltyp)
835 CASE DEFAULT
836 cpabort("")
837 CASE (do_nonrel_atom, do_zoramp_atom, do_sczoramp_atom)
838 dkhorder = -1
839 CASE (do_dkh0_atom)
840 dkhorder = 0
841 CASE (do_dkh1_atom)
842 dkhorder = 1
843 CASE (do_dkh2_atom)
844 dkhorder = 2
845 CASE (do_dkh3_atom)
846 dkhorder = 3
847 END SELECT
848
849 SELECT CASE (reltyp)
850 CASE DEFAULT
851 cpabort("")
852 CASE (do_nonrel_atom)
853 ! nothing to do
854 NULLIFY (integrals%tzora, integrals%hdkh)
855 CASE (do_zoramp_atom, do_sczoramp_atom)
856
857 NULLIFY (integrals%hdkh)
858
859 IF (integrals%zorastat == 0) THEN
860 n = maxval(basis%nbas)
861 ALLOCATE (integrals%tzora(n, n, 0:lmat))
862 integrals%tzora = 0._dp
863 m = basis%grid%nr
864 ALLOCATE (modpot(1:m), cpot(1:m))
865 CALL calculate_model_potential(modpot, basis%grid, zcore)
866 ! Zora potential
867 cpot(1:m) = modpot(1:m)/(4._dp*c_light_au*c_light_au - 2._dp*modpot(1:m))
868 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
869 CALL numpot_matrix(integrals%tzora, cpot, basis, 0)
870 DO l = 0, lmat
871 nl = basis%nbas(l)
872 integrals%tzora(1:nl, 1:nl, l) = real(l*(l + 1), dp)*integrals%tzora(1:nl, 1:nl, l)
873 END DO
874 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
875 CALL numpot_matrix(integrals%tzora, cpot, basis, 2)
876 !
877 ! scaled ZORA
878 IF (reltyp == do_sczoramp_atom) THEN
879 ALLOCATE (hmat(n, n, 0:lmat), wfn(n, n, 0:lmat), ener(n, 0:lmat), pvp(n, n, 0:lmat), sps(n, n))
880 hmat(:, :, :) = integrals%kin + integrals%tzora
881 ! model potential
882 CALL numpot_matrix(hmat, modpot, basis, 0)
883 ! eigenvalues and eigenvectors
884 CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne, lmat)
885 ! relativistic kinetic energy
886 cpot(1:m) = c_light_au*c_light_au/(2._dp*c_light_au*c_light_au - modpot(1:m))**2
887 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
888 pvp = 0.0_dp
889 CALL numpot_matrix(pvp, cpot, basis, 0)
890 DO l = 0, lmat
891 nl = basis%nbas(l)
892 pvp(1:nl, 1:nl, l) = real(l*(l + 1), dp)*pvp(1:nl, 1:nl, l)
893 END DO
894 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
895 CALL numpot_matrix(pvp, cpot, basis, 2)
896 ! calculate psi*pvp*psi and the scaled orbital energies
897 ! actually, we directly calculate the energy difference
898 DO l = 0, lmat
899 nl = basis%nbas(l)
900 DO i = 1, integrals%nne(l)
901 IF (ener(i, l) < 0._dp) THEN
902 ascal = sum(wfn(1:nl, i, l)*matmul(pvp(1:nl, 1:nl, l), wfn(1:nl, i, l)))
903 ener(i, l) = ener(i, l)*ascal/(1.0_dp + ascal)
904 ELSE
905 ener(i, l) = 0.0_dp
906 END IF
907 END DO
908 END DO
909 ! correction term is calculated as a projector
910 hmat = 0.0_dp
911 DO l = 0, lmat
912 nl = basis%nbas(l)
913 DO i = 1, integrals%nne(l)
914 DO k1 = 1, nl
915 DO k2 = 1, nl
916 hmat(k1, k2, l) = hmat(k1, k2, l) + ener(i, l)*wfn(k1, i, l)*wfn(k2, i, l)
917 END DO
918 END DO
919 END DO
920 ! transform with overlap matrix
921 sps(1:nl, 1:nl) = matmul(integrals%ovlp(1:nl, 1:nl, l), &
922 matmul(hmat(1:nl, 1:nl, l), integrals%ovlp(1:nl, 1:nl, l)))
923 ! add scaling correction to tzora
924 integrals%tzora(1:nl, 1:nl, l) = integrals%tzora(1:nl, 1:nl, l) - sps(1:nl, 1:nl)
925 END DO
926
927 DEALLOCATE (hmat, wfn, ener, pvp, sps)
928 END IF
929 !
930 DEALLOCATE (modpot, cpot)
931
932 integrals%zorastat = 1
933
934 END IF
935
936 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
937
938 NULLIFY (integrals%tzora)
939
940 IF (integrals%dkhstat == 0) THEN
941 n = maxval(basis%nbas)
942 ALLOCATE (integrals%hdkh(n, n, 0:lmat))
943 integrals%hdkh = 0._dp
944
945 m = maxval(basis%nprim)
946 ALLOCATE (tp(m, m, 0:lmat), sp(m, m, 0:lmat), vp(m, m, 0:lmat), pvp(m, m, 0:lmat))
947 tp = 0._dp
948 sp = 0._dp
949 vp = 0._dp
950 pvp = 0._dp
951
952 SELECT CASE (basis%basis_type)
953 CASE DEFAULT
954 cpabort("")
955 CASE (gto_basis, cgto_basis)
956
957 DO l = 0, lmat
958 m = basis%nprim(l)
959 IF (m > 0) THEN
960 CALL sg_kinetic(tp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
961 CALL sg_overlap(sp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
962 IF (PRESENT(alpha)) THEN
963 CALL sg_erfc(vp(1:m, 1:m, l), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
964 ELSE
965 CALL sg_nuclear(vp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
966 END IF
967 CALL sg_kinnuc(pvp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
968 vp(1:m, 1:m, l) = -zcore*vp(1:m, 1:m, l)
969 pvp(1:m, 1:m, l) = -zcore*pvp(1:m, 1:m, l)
970 END IF
971 END DO
972
973 CASE (sto_basis)
974 cpabort("")
975 CASE (num_basis)
976 cpabort("")
977 END SELECT
978
979 CALL dkh_integrals(integrals, basis, dkhorder, sp, tp, vp, pvp)
980
981 integrals%dkhstat = 1
982
983 DEALLOCATE (tp, sp, vp, pvp)
984
985 ELSE
986 cpassert(ASSOCIATED(integrals%hdkh))
987 END IF
988
989 END SELECT
990
991 CALL timestop(handle)
992
993 END SUBROUTINE atom_relint_setup
994
995! **************************************************************************************************
996!> \brief ...
997!> \param integrals ...
998!> \param basis ...
999!> \param order ...
1000!> \param sp ...
1001!> \param tp ...
1002!> \param vp ...
1003!> \param pvp ...
1004! **************************************************************************************************
1005 SUBROUTINE dkh_integrals(integrals, basis, order, sp, tp, vp, pvp)
1006 TYPE(atom_integrals), INTENT(INOUT) :: integrals
1007 TYPE(atom_basis_type), INTENT(INOUT) :: basis
1008 INTEGER, INTENT(IN) :: order
1009 REAL(dp), DIMENSION(:, :, 0:) :: sp, tp, vp, pvp
1010
1011 INTEGER :: l, m, n
1012 REAL(dp), DIMENSION(:, :, :), POINTER :: hdkh
1013
1014 cpassert(order >= 0)
1015
1016 hdkh => integrals%hdkh
1017
1018 DO l = 0, lmat
1019 n = integrals%n(l)
1020 m = basis%nprim(l)
1021 IF (n > 0) THEN
1022 CALL dkh_atom_transformation(sp(1:m, 1:m, l), vp(1:m, 1:m, l), tp(1:m, 1:m, l), pvp(1:m, 1:m, l), m, order)
1023 SELECT CASE (basis%basis_type)
1024 CASE DEFAULT
1025 cpabort("")
1026 CASE (gto_basis)
1027 cpassert(n == m)
1028 integrals%hdkh(1:n, 1:n, l) = tp(1:n, 1:n, l) + vp(1:n, 1:n, l)
1029 CASE (cgto_basis)
1030 CALL contract2(integrals%hdkh(1:n, 1:n, l), tp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1031 CALL contract2add(integrals%hdkh(1:n, 1:n, l), vp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1032 CASE (sto_basis)
1033 cpabort("")
1034 CASE (num_basis)
1035 cpabort("")
1036 END SELECT
1037 ELSE
1038 integrals%hdkh(1:n, 1:n, l) = 0._dp
1039 END IF
1040 END DO
1041
1042 END SUBROUTINE dkh_integrals
1043
1044! **************************************************************************************************
1045!> \brief Calculate overlap matrix between two atomic basis sets.
1046!> \param ovlap overlap matrix <chi_{a,l} | chi_{b,l}>
1047!> \param basis_a first basis set
1048!> \param basis_b second basis set
1049! **************************************************************************************************
1050 SUBROUTINE atom_basis_projection_overlap(ovlap, basis_a, basis_b)
1051 REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: ovlap
1052 TYPE(atom_basis_type), INTENT(IN) :: basis_a, basis_b
1053
1054 INTEGER :: i, j, l, ma, mb, na, nb
1055 LOGICAL :: ebas
1056 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
1057
1058 ebas = (basis_a%basis_type == basis_b%basis_type)
1059
1060 ovlap = 0.0_dp
1061
1062 IF (ebas) THEN
1063 SELECT CASE (basis_a%basis_type)
1064 CASE DEFAULT
1065 cpabort("")
1066 CASE (gto_basis)
1067 DO l = 0, lmat
1068 na = basis_a%nbas(l)
1069 nb = basis_b%nbas(l)
1070 CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l))
1071 END DO
1072 CASE (cgto_basis)
1073 DO l = 0, lmat
1074 na = basis_a%nbas(l)
1075 nb = basis_b%nbas(l)
1076 ma = basis_a%nprim(l)
1077 mb = basis_b%nprim(l)
1078 ALLOCATE (omat(ma, mb))
1079 CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l))
1080 ovlap(1:na, 1:nb, l) = matmul(transpose(basis_a%cm(1:ma, 1:na, l)), &
1081 matmul(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l)))
1082 DEALLOCATE (omat)
1083 END DO
1084 CASE (sto_basis)
1085 DO l = 0, lmat
1086 na = basis_a%nbas(l)
1087 nb = basis_b%nbas(l)
1088 CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), &
1089 basis_a%ns(1:na, l), basis_b%as(1:nb, l))
1090 END DO
1091 CASE (num_basis)
1092 cpassert(atom_compare_grids(basis_a%grid, basis_b%grid))
1093 DO l = 0, lmat
1094 na = basis_a%nbas(l)
1095 nb = basis_b%nbas(l)
1096 DO j = 1, nb
1097 DO i = 1, na
1098 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1099 END DO
1100 END DO
1101 END DO
1102 END SELECT
1103 ELSE
1104 cpassert(atom_compare_grids(basis_a%grid, basis_b%grid))
1105 DO l = 0, lmat
1106 na = basis_a%nbas(l)
1107 nb = basis_b%nbas(l)
1108 DO j = 1, nb
1109 DO i = 1, na
1110 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1111 END DO
1112 END DO
1113 END DO
1114 END IF
1115
1116 END SUBROUTINE atom_basis_projection_overlap
1117
1118! **************************************************************************************************
1119!> \brief Release memory allocated for atomic integrals (valence electrons).
1120!> \param integrals atomic integrals
1121! **************************************************************************************************
1122 SUBROUTINE atom_int_release(integrals)
1123 TYPE(atom_integrals), INTENT(INOUT) :: integrals
1124
1125 INTEGER :: ll
1126
1127 IF (ASSOCIATED(integrals%ovlp)) THEN
1128 DEALLOCATE (integrals%ovlp)
1129 END IF
1130 IF (ASSOCIATED(integrals%kin)) THEN
1131 DEALLOCATE (integrals%kin)
1132 END IF
1133 IF (ASSOCIATED(integrals%conf)) THEN
1134 DEALLOCATE (integrals%conf)
1135 END IF
1136 DO ll = 1, SIZE(integrals%ceri)
1137 IF (ASSOCIATED(integrals%ceri(ll)%int)) THEN
1138 DEALLOCATE (integrals%ceri(ll)%int)
1139 END IF
1140 IF (ASSOCIATED(integrals%eeri(ll)%int)) THEN
1141 DEALLOCATE (integrals%eeri(ll)%int)
1142 END IF
1143 END DO
1144 IF (ASSOCIATED(integrals%utrans)) THEN
1145 DEALLOCATE (integrals%utrans)
1146 END IF
1147 IF (ASSOCIATED(integrals%uptrans)) THEN
1148 DEALLOCATE (integrals%uptrans)
1149 END IF
1150
1151 integrals%status = 0
1152
1153 END SUBROUTINE atom_int_release
1154
1155! **************************************************************************************************
1156!> \brief Release memory allocated for atomic integrals (core electrons).
1157!> \param integrals atomic integrals
1158! **************************************************************************************************
1159 SUBROUTINE atom_ppint_release(integrals)
1160 TYPE(atom_integrals), INTENT(INOUT) :: integrals
1161
1162 IF (ASSOCIATED(integrals%hnl)) THEN
1163 DEALLOCATE (integrals%hnl)
1164 END IF
1165 IF (ASSOCIATED(integrals%core)) THEN
1166 DEALLOCATE (integrals%core)
1167 END IF
1168 IF (ASSOCIATED(integrals%clsd)) THEN
1169 DEALLOCATE (integrals%clsd)
1170 END IF
1171
1172 integrals%ppstat = 0
1173
1174 END SUBROUTINE atom_ppint_release
1175
1176! **************************************************************************************************
1177!> \brief Release memory allocated for atomic integrals (relativistic effects).
1178!> \param integrals atomic integrals
1179! **************************************************************************************************
1180 SUBROUTINE atom_relint_release(integrals)
1181 TYPE(atom_integrals), INTENT(INOUT) :: integrals
1182
1183 IF (ASSOCIATED(integrals%tzora)) THEN
1184 DEALLOCATE (integrals%tzora)
1185 END IF
1186 IF (ASSOCIATED(integrals%hdkh)) THEN
1187 DEALLOCATE (integrals%hdkh)
1188 END IF
1189
1190 integrals%zorastat = 0
1191 integrals%dkhstat = 0
1192
1193 END SUBROUTINE atom_relint_release
1194
1195! **************************************************************************************************
1196!> \brief Calculate model potential. It is useful to guess atomic orbitals.
1197!> \param modpot model potential
1198!> \param grid atomic radial grid
1199!> \param zcore nuclear charge
1200! **************************************************************************************************
1201 SUBROUTINE calculate_model_potential(modpot, grid, zcore)
1202 REAL(dp), DIMENSION(:), INTENT(INOUT) :: modpot
1203 TYPE(grid_atom_type), INTENT(IN) :: grid
1204 REAL(dp), INTENT(IN) :: zcore
1205
1206 INTEGER :: i, ii, l, ll, n, z
1207 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: pot, rho
1208 TYPE(atom_state) :: state
1209
1210 n = SIZE(modpot)
1211 ALLOCATE (rho(n), pot(n))
1212
1213 ! fill default occupation
1214 state%core = 0._dp
1215 state%occ = 0._dp
1216 state%multiplicity = -1
1217 z = nint(zcore)
1218 DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
1219 IF (ptable(z)%e_conv(l) /= 0) THEN
1220 state%maxl_occ = l
1221 ll = 2*(2*l + 1)
1222 DO i = 1, SIZE(state%occ, 2)
1223 ii = ptable(z)%e_conv(l) - (i - 1)*ll
1224 IF (ii <= ll) THEN
1225 state%occ(l, i) = ii
1226 EXIT
1227 ELSE
1228 state%occ(l, i) = ll
1229 END IF
1230 END DO
1231 END IF
1232 END DO
1233
1234 modpot = -zcore/grid%rad(:)
1235
1236 ! Coulomb potential
1237 CALL slater_density(rho, pot, nint(zcore), state, grid)
1238 CALL coulomb_potential_numeric(pot, rho, grid)
1239 modpot = modpot + pot
1240
1241 ! XC potential
1242 CALL wigner_slater_functional(rho, pot)
1243 modpot = modpot + pot
1244
1245 DEALLOCATE (rho, pot)
1246
1247 END SUBROUTINE calculate_model_potential
1248
1249END MODULE atom_operators
atom_utils::integrate_grid
Definition atom_utils.F:82
ai_onecenter
Definition ai_onecenter.F:28
ai_onecenter::sg_kinnuc
subroutine, public sg_kinnuc(umat, l, pa, pb)
...
Definition ai_onecenter.F:172
ai_onecenter::sg_nuclear
subroutine, public sg_nuclear(umat, l, pa, pb)
...
Definition ai_onecenter.F:137
ai_onecenter::sg_coulomb
subroutine, public sg_coulomb(eri, nu, pa, lab, pc, lcd)
...
Definition ai_onecenter.F:452
ai_onecenter::sg_exchange
subroutine, public sg_exchange(eri, nu, pa, lac, pb, lbd)
...
Definition ai_onecenter.F:516
ai_onecenter::sg_erfc
subroutine, public sg_erfc(umat, l, a, pa, pb)
...
Definition ai_onecenter.F:590
ai_onecenter::sg_kinetic
subroutine, public sg_kinetic(kmat, l, pa, pb)
...
Definition ai_onecenter.F:102
ai_onecenter::sto_nuclear
subroutine, public sto_nuclear(umat, na, pa, nb, pb)
...
Definition ai_onecenter.F:757
ai_onecenter::sto_overlap
subroutine, public sto_overlap(smat, na, pa, nb, pb)
...
Definition ai_onecenter.F:668
ai_onecenter::sto_kinetic
subroutine, public sto_kinetic(kmat, l, na, pa, nb, pb)
...
Definition ai_onecenter.F:711
ai_onecenter::sg_gpot
subroutine, public sg_gpot(vpmat, k, rc, l, pa, pb)
...
Definition ai_onecenter.F:373
ai_onecenter::sg_proj_ol
subroutine, public sg_proj_ol(spmat, l, p, k, rc)
...
Definition ai_onecenter.F:337
ai_onecenter::sg_erf
subroutine, public sg_erf(upmat, l, a, pa, pb)
...
Definition ai_onecenter.F:231
ai_onecenter::sg_overlap
subroutine, public sg_overlap(smat, l, pa, pb)
...
Definition ai_onecenter.F:65
atom_operators
Calculate the atomic operator matrices.
Definition atom_operators.F:15
atom_operators::atom_ppint_release
subroutine, public atom_ppint_release(integrals)
Release memory allocated for atomic integrals (core electrons).
Definition atom_operators.F:1160
atom_operators::atom_int_setup
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
Definition atom_operators.F:63
atom_operators::atom_relint_setup
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
Definition atom_operators.F:817
atom_operators::atom_relint_release
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
Definition atom_operators.F:1181
atom_operators::atom_basis_projection_overlap
subroutine, public atom_basis_projection_overlap(ovlap, basis_a, basis_b)
Calculate overlap matrix between two atomic basis sets.
Definition atom_operators.F:1051
atom_operators::atom_ppint_setup
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
Definition atom_operators.F:336
atom_operators::calculate_model_potential
subroutine, public calculate_model_potential(modpot, grid, zcore)
Calculate model potential. It is useful to guess atomic orbitals.
Definition atom_operators.F:1202
atom_operators::atom_int_release
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Definition atom_operators.F:1123
atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
atom_types::num_basis
integer, parameter, public num_basis
Definition atom_types.F:69
atom_types::cgto_basis
integer, parameter, public cgto_basis
Definition atom_types.F:69
atom_types::gto_basis
integer, parameter, public gto_basis
Definition atom_types.F:69
atom_types::sto_basis
integer, parameter, public sto_basis
Definition atom_types.F:69
atom_types::atom_compare_grids
logical function, public atom_compare_grids(grid1, grid2)
...
Definition atom_types.F:3065
atom_types::lmat
integer, parameter, public lmat
Definition atom_types.F:67
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:913
atom_types::atom_basis_gridrep
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
Definition atom_types.F:781
atom_utils
Some basic routines for atomic calculations.
Definition atom_utils.F:15
atom_utils::slater_density
subroutine, public slater_density(density1, density2, zcore, state, grid)
Calculate Slater density on a radial grid.
Definition atom_utils.F:2093
atom_utils::contract2
subroutine, public contract2(int, omat, cm)
Transform a matrix expressed in terms of a uncontracted basis set to a contracted one.
Definition atom_utils.F:2612
atom_utils::wigner_slater_functional
pure subroutine, public wigner_slater_functional(rho, vxc)
Calculate the functional derivative of the Wigner (correlation) - Slater (exchange) density functiona...
Definition atom_utils.F:2158
atom_utils::contract4
subroutine, public contract4(eri, omat, cm1, cm2)
Contract a matrix of Electron Repulsion Integrals (ERI-s).
Definition atom_utils.F:2664
atom_utils::numpot_matrix
subroutine, public numpot_matrix(imat, cpot, basis, derivatives)
Calculate a potential matrix by integrating the potential on an atomic radial grid.
Definition atom_utils.F:1839
atom_utils::atom_solve
subroutine, public atom_solve(hmat, umat, orb, ener, nb, nv, maxl)
Solve the generalised eigenproblem for every angular momentum.
Definition atom_utils.F:943
atom_utils::contract2add
subroutine, public contract2add(int, omat, cm)
Same as contract2(), but add the new contracted matrix to the old one instead of overwriting it.
Definition atom_utils.F:2638
atom_utils::coulomb_potential_numeric
subroutine, public coulomb_potential_numeric(cpot, density, grid)
Numerically compute the Coulomb potential on an atomic radial grid.
Definition atom_utils.F:1076
dkh_main
Definition dkh_main.F:7
dkh_main::dkh_atom_transformation
subroutine, public dkh_atom_transformation(s, v, h, pvp, n, dkh_order)
...
Definition dkh_main.F:1245
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::sgp_pseudo
integer, parameter, public sgp_pseudo
Definition input_constants.F:174
input_constants::do_dkh3_atom
integer, parameter, public do_dkh3_atom
Definition input_constants.F:165
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:174
input_constants::ecp_pseudo
integer, parameter, public ecp_pseudo
Definition input_constants.F:174
input_constants::do_nonrel_atom
integer, parameter, public do_nonrel_atom
Definition input_constants.F:165
input_constants::do_dkh0_atom
integer, parameter, public do_dkh0_atom
Definition input_constants.F:165
input_constants::upf_pseudo
integer, parameter, public upf_pseudo
Definition input_constants.F:174
input_constants::poly_conf
integer, parameter, public poly_conf
Definition input_constants.F:180
input_constants::do_dkh2_atom
integer, parameter, public do_dkh2_atom
Definition input_constants.F:165
input_constants::no_pseudo
integer, parameter, public no_pseudo
Definition input_constants.F:174
input_constants::barrier_conf
integer, parameter, public barrier_conf
Definition input_constants.F:180
input_constants::do_zoramp_atom
integer, parameter, public do_zoramp_atom
Definition input_constants.F:165
input_constants::do_dkh1_atom
integer, parameter, public do_dkh1_atom
Definition input_constants.F:165
input_constants::do_sczoramp_atom
integer, parameter, public do_sczoramp_atom
Definition input_constants.F:165
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::gamma1
real(kind=dp), dimension(0:maxfac), parameter, public gamma1
Definition mathconstants.F:95
mathconstants::sqrt2
real(kind=dp), parameter, public sqrt2
Definition mathconstants.F:127
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::jacobi
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition mathlib.F:1468
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::c_light_au
real(kind=dp), parameter, public c_light_au
Definition physcon.F:90
qs_grid_atom
Definition qs_grid_atom.F:8
atom_types::atom_basis_type
Provides all information about a basis set.
Definition atom_types.F:78
atom_types::atom_integrals
Definition atom_types.F:216
atom_types::atom_potential_type
Definition atom_types.F:171
atom_types::atom_state
Provides all information on states and occupation.
Definition atom_types.F:198
atom_types::eri
Holds atomic integrals.
Definition atom_types.F:212
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44