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