(git:b77b4be)
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 alpha = 1._dp/potential%gth_pot%rc/sqrt(2._dp)
380 DO l = 0, lmat
381 n = integrals%n(l)
382 ALLOCATE (omat(n, n), spmat(n, 5))
383
384 omat = 0._dp
385 CALL sg_erf(omat(1:n, 1:n), l, alpha, basis%am(1:n, l), basis%am(1:n, l))
386 integrals%core(1:n, 1:n, l) = -potential%gth_pot%zion*omat(1:n, 1:n)
387 DO i = 1, potential%gth_pot%ncl
388 omat = 0._dp
389 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))
390 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
391 potential%gth_pot%cl(i)*omat(1:n, 1:n)
392 END DO
393 IF (potential%gth_pot%lpotextended) THEN
394 DO k = 1, potential%gth_pot%nexp_lpot
395 DO i = 1, potential%gth_pot%nct_lpot(k)
396 omat = 0._dp
397 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lpot(k), l, &
398 basis%am(1:n, l), basis%am(1:n, l))
399 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
400 potential%gth_pot%cval_lpot(i, k)*omat(1:n, 1:n)
401 END DO
402 END DO
403 END IF
404 IF (potential%gth_pot%lsdpot) THEN
405 DO k = 1, potential%gth_pot%nexp_lsd
406 DO i = 1, potential%gth_pot%nct_lsd(k)
407 omat = 0._dp
408 CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lsd(k), l, &
409 basis%am(1:n, l), basis%am(1:n, l))
410 integrals%clsd(1:n, 1:n, l) = integrals%clsd(1:n, 1:n, l) + &
411 potential%gth_pot%cval_lsd(i, k)*omat(1:n, 1:n)
412 END DO
413 END DO
414 END IF
415
416 spmat = 0._dp
417 m = potential%gth_pot%nl(l)
418 DO i = 1, m
419 CALL sg_proj_ol(spmat(1:n, i), l, basis%am(1:n, l), i - 1, potential%gth_pot%rcnl(l))
420 END DO
421 integrals%hnl(1:n, 1:n, l) = matmul(spmat(1:n, 1:m), &
422 matmul(potential%gth_pot%hnl(1:m, 1:m, l), transpose(spmat(1:n, 1:m))))
423
424 DEALLOCATE (omat, spmat)
425 END DO
426 CASE (upf_pseudo)
427 CALL upfint_setup(integrals, basis, potential)
428 CASE (sgp_pseudo)
429 CALL sgpint_setup(integrals, basis, potential)
430 CASE DEFAULT
431 cpabort("")
432 END SELECT
433
434 CASE (cgto_basis)
435
436 SELECT CASE (potential%ppot_type)
437 CASE (no_pseudo, ecp_pseudo)
438 DO l = 0, lmat
439 n = integrals%n(l)
440 m = basis%nprim(l)
441 ALLOCATE (omat(m, m))
442
443 CALL sg_nuclear(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
444 CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
445
446 DEALLOCATE (omat)
447 END DO
448 CASE (gth_pseudo)
449 alpha = 1._dp/potential%gth_pot%rc/sqrt(2._dp)
450 DO l = 0, lmat
451 n = integrals%n(l)
452 m = basis%nprim(l)
453 IF (n > 0 .AND. m > 0) THEN
454 ALLOCATE (omat(m, m), spmat(n, 5), xmat(m))
455
456 omat = 0._dp
457 CALL sg_erf(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
458 omat(1:m, 1:m) = -potential%gth_pot%zion*omat(1:m, 1:m)
459 CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
460 DO i = 1, potential%gth_pot%ncl
461 omat = 0._dp
462 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))
463 omat(1:m, 1:m) = potential%gth_pot%cl(i)*omat(1:m, 1:m)
464 CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
465 END DO
466 IF (potential%gth_pot%lpotextended) THEN
467 DO k = 1, potential%gth_pot%nexp_lpot
468 DO i = 1, potential%gth_pot%nct_lpot(k)
469 omat = 0._dp
470 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lpot(k), l, &
471 basis%am(1:m, l), basis%am(1:m, l))
472 omat(1:m, 1:m) = potential%gth_pot%cval_lpot(i, k)*omat(1:m, 1:m)
473 CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
474 END DO
475 END DO
476 END IF
477 IF (potential%gth_pot%lsdpot) THEN
478 DO k = 1, potential%gth_pot%nexp_lsd
479 DO i = 1, potential%gth_pot%nct_lsd(k)
480 omat = 0._dp
481 CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lsd(k), l, &
482 basis%am(1:m, l), basis%am(1:m, l))
483 omat(1:m, 1:m) = potential%gth_pot%cval_lsd(i, k)*omat(1:m, 1:m)
484 CALL contract2add(integrals%clsd(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
485 END DO
486 END DO
487 END IF
488
489 spmat = 0._dp
490 k = potential%gth_pot%nl(l)
491 DO i = 1, k
492 CALL sg_proj_ol(xmat(1:m), l, basis%am(1:m, l), i - 1, potential%gth_pot%rcnl(l))
493 spmat(1:n, i) = matmul(transpose(basis%cm(1:m, 1:n, l)), xmat(1:m))
494 END DO
495 IF (k > 0) THEN
496 integrals%hnl(1:n, 1:n, l) = matmul(spmat(1:n, 1:k), &
497 matmul(potential%gth_pot%hnl(1:k, 1:k, l), &
498 transpose(spmat(1:n, 1:k))))
499 END IF
500 DEALLOCATE (omat, spmat, xmat)
501 END IF
502 END DO
503 CASE (upf_pseudo)
504 CALL upfint_setup(integrals, basis, potential)
505 CASE (sgp_pseudo)
506 CALL sgpint_setup(integrals, basis, potential)
507 CASE DEFAULT
508 cpabort("")
509 END SELECT
510
511 CASE (sto_basis)
512
513 SELECT CASE (potential%ppot_type)
514 CASE (no_pseudo, ecp_pseudo)
515 DO l = 0, lmat
516 n = integrals%n(l)
517 CALL sto_nuclear(integrals%core(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
518 basis%ns(1:n, l), basis%as(1:n, l))
519 END DO
520 CASE (gth_pseudo)
521 rad => basis%grid%rad
522 m = basis%grid%nr
523 ALLOCATE (cpot(1:m))
524 rc = potential%gth_pot%rc
525 alpha = 1._dp/rc/sqrt(2._dp)
526
527 ! local pseudopotential, we use erf = 1/r - erfc
528 integrals%core = 0._dp
529 DO i = 1, m
530 cpot(i) = potential%gth_pot%zion*erfc(alpha*rad(i))/rad(i)
531 END DO
532 DO i = 1, potential%gth_pot%ncl
533 ii = 2*(i - 1)
534 cpot(1:m) = cpot(1:m) + potential%gth_pot%cl(i)*(rad/rc)**ii*exp(-0.5_dp*(rad/rc)**2)
535 END DO
536 IF (potential%gth_pot%lpotextended) THEN
537 DO k = 1, potential%gth_pot%nexp_lpot
538 al = potential%gth_pot%alpha_lpot(k)
539 DO i = 1, potential%gth_pot%nct_lpot(k)
540 ii = 2*(i - 1)
541 cpot(1:m) = cpot(1:m) + potential%gth_pot%cval_lpot(i, k)*(rad/al)**ii*exp(-0.5_dp*(rad/al)**2)
542 END DO
543 END DO
544 END IF
545 CALL numpot_matrix(integrals%core, cpot, basis, 0)
546 DO l = 0, lmat
547 n = integrals%n(l)
548 ALLOCATE (omat(n, n))
549 omat = 0._dp
550 CALL sto_nuclear(omat(1:n, 1:n), basis%ns(1:n, l), basis%as(1:n, l), &
551 basis%ns(1:n, l), basis%as(1:n, l))
552 integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) - potential%gth_pot%zion*omat(1:n, 1:n)
553 DEALLOCATE (omat)
554 END DO
555
556 IF (potential%gth_pot%lsdpot) THEN
557 cpot = 0._dp
558 DO k = 1, potential%gth_pot%nexp_lsd
559 al = potential%gth_pot%alpha_lsd(k)
560 DO i = 1, potential%gth_pot%nct_lsd(k)
561 ii = 2*(i - 1)
562 cpot(:) = cpot + potential%gth_pot%cval_lsd(i, k)*(rad/al)**ii*exp(-0.5_dp*(rad/al)**2)
563 END DO
564 END DO
565 CALL numpot_matrix(integrals%clsd, cpot, basis, 0)
566 END IF
567
568 DO l = 0, lmat
569 n = integrals%n(l)
570 ! non local pseudopotential
571 ALLOCATE (spmat(n, 5))
572 spmat = 0._dp
573 k = potential%gth_pot%nl(l)
574 DO i = 1, k
575 rc = potential%gth_pot%rcnl(l)
576 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)
577 DO j = 1, basis%nbas(l)
578 spmat(j, i) = integrate_grid(cpot, basis%bf(:, j, l), basis%grid)
579 END DO
580 END DO
581 integrals%hnl(1:n, 1:n, l) = matmul(spmat(1:n, 1:k), &
582 matmul(potential%gth_pot%hnl(1:k, 1:k, l), &
583 transpose(spmat(1:n, 1:k))))
584 DEALLOCATE (spmat)
585 END DO
586 DEALLOCATE (cpot)
587
588 CASE (upf_pseudo)
589 CALL upfint_setup(integrals, basis, potential)
590 CASE (sgp_pseudo)
591 CALL sgpint_setup(integrals, basis, potential)
592 CASE DEFAULT
593 cpabort("")
594 END SELECT
595
596 CASE (num_basis)
597 cpabort("")
598 END SELECT
599
600 ! add ecp_pseudo using numerical representation of basis
601 IF (potential%ppot_type == ecp_pseudo) THEN
602 ! scale 1/r potential
603 integrals%core = -potential%ecp_pot%zion*integrals%core
604 ! local potential
605 m = basis%grid%nr
606 rad => basis%grid%rad
607 ALLOCATE (cpot(1:m))
608 cpot = 0._dp
609 DO k = 1, potential%ecp_pot%nloc
610 n = potential%ecp_pot%nrloc(k)
611 alpha = potential%ecp_pot%bloc(k)
612 cpot(:) = cpot + potential%ecp_pot%aloc(k)*rad**(n - 2)*exp(-alpha*rad**2)
613 END DO
614 CALL numpot_matrix(integrals%core, cpot, basis, 0)
615 ! non local pseudopotential
616 DO l = 0, min(potential%ecp_pot%lmax, lmat)
617 cpot = 0._dp
618 DO k = 1, potential%ecp_pot%npot(l)
619 n = potential%ecp_pot%nrpot(k, l)
620 alpha = potential%ecp_pot%bpot(k, l)
621 cpot(:) = cpot + potential%ecp_pot%apot(k, l)*rad**(n - 2)*exp(-alpha*rad**2)
622 END DO
623 DO i = 1, basis%nbas(l)
624 DO j = i, basis%nbas(l)
625 integrals%hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
626 integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
627 END DO
628 END DO
629 END DO
630 DEALLOCATE (cpot)
631 END IF
632
633 END IF
634
635 CALL timestop(handle)
636
637 END SUBROUTINE atom_ppint_setup
638
639! **************************************************************************************************
640!> \brief ...
641!> \param integrals ...
642!> \param basis ...
643!> \param potential ...
644! **************************************************************************************************
645 SUBROUTINE upfint_setup(integrals, basis, potential)
646 TYPE(atom_integrals), INTENT(INOUT) :: integrals
647 TYPE(atom_basis_type), INTENT(INOUT) :: basis
648 TYPE(atom_potential_type), INTENT(IN) :: potential
649
650 CHARACTER(len=4) :: ptype
651 INTEGER :: i, j, k1, k2, la, lb, m, n
652 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: spot
653 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: spmat
654 TYPE(atom_basis_type) :: gbasis
655
656 ! get basis representation on UPF grid
657 CALL atom_basis_gridrep(basis, gbasis, potential%upf_pot%r, potential%upf_pot%rab)
658
659 ! local pseudopotential
660 integrals%core = 0._dp
661 CALL numpot_matrix(integrals%core, potential%upf_pot%vlocal, gbasis, 0)
662
663 ptype = adjustl(trim(potential%upf_pot%pseudo_type))
664 IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
665 ! non local pseudopotential
666 n = maxval(integrals%n(:))
667 m = potential%upf_pot%number_of_proj
668 ALLOCATE (spmat(n, m))
669 spmat = 0.0_dp
670 DO i = 1, m
671 la = potential%upf_pot%lbeta(i)
672 DO j = 1, gbasis%nbas(la)
673 spmat(j, i) = integrate_grid(potential%upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
674 END DO
675 END DO
676 DO i = 1, m
677 la = potential%upf_pot%lbeta(i)
678 DO j = 1, m
679 lb = potential%upf_pot%lbeta(j)
680 IF (la == lb) THEN
681 DO k1 = 1, gbasis%nbas(la)
682 DO k2 = 1, gbasis%nbas(la)
683 integrals%hnl(k1, k2, la) = integrals%hnl(k1, k2, la) + &
684 spmat(k1, i)*potential%upf_pot%dion(i, j)*spmat(k2, j)
685 END DO
686 END DO
687 END IF
688 END DO
689 END DO
690 DEALLOCATE (spmat)
691 ELSE IF (ptype(1:2) == "SL") THEN
692 ! semi local pseudopotential
693 DO la = 0, potential%upf_pot%l_max
694 IF (la == potential%upf_pot%l_local) cycle
695 m = SIZE(potential%upf_pot%vsemi(:, la + 1))
696 ALLOCATE (spot(m))
697 spot(:) = potential%upf_pot%vsemi(:, la + 1) - potential%upf_pot%vlocal(:)
698 n = basis%nbas(la)
699 DO i = 1, n
700 DO j = i, n
701 integrals%core(i, j, la) = integrals%core(i, j, la) + &
702 integrate_grid(spot(:), &
703 gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
704 integrals%core(j, i, la) = integrals%core(i, j, la)
705 END DO
706 END DO
707 DEALLOCATE (spot)
708 END DO
709 ELSE
710 cpabort("Pseudopotential type: ["//adjustl(trim(ptype))//"] not known")
711 END IF
712
713 ! release basis representation on UPF grid
714 CALL release_atom_basis(gbasis)
715
716 END SUBROUTINE upfint_setup
717
718! **************************************************************************************************
719!> \brief ...
720!> \param integrals ...
721!> \param basis ...
722!> \param potential ...
723! **************************************************************************************************
724 SUBROUTINE sgpint_setup(integrals, basis, potential)
725 TYPE(atom_integrals), INTENT(INOUT) :: integrals
726 TYPE(atom_basis_type), INTENT(INOUT) :: basis
727 TYPE(atom_potential_type), INTENT(IN) :: potential
728
729 INTEGER :: i, ia, j, l, m, n, na
730 REAL(kind=dp) :: a, c, rc, zval
731 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot, pgauss
732 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: qmat
733 REAL(kind=dp), DIMENSION(:), POINTER :: rad
734
735 rad => basis%grid%rad
736 m = basis%grid%nr
737
738 ! local pseudopotential
739 integrals%core = 0._dp
740 ALLOCATE (cpot(m))
741 cpot = 0.0_dp
742 zval = potential%sgp_pot%zion
743 DO i = 1, m
744 rc = rad(i)/potential%sgp_pot%ac_local/sqrt(2.0_dp)
745 cpot(i) = cpot(i) - zval/rad(i)*erf(rc)
746 END DO
747 DO i = 1, potential%sgp_pot%n_local
748 cpot(:) = cpot(:) + potential%sgp_pot%c_local(i)*exp(-potential%sgp_pot%a_local(i)*rad(:)**2)
749 END DO
750 CALL numpot_matrix(integrals%core, cpot, basis, 0)
751 DEALLOCATE (cpot)
752
753 ! nonlocal pseudopotential
754 integrals%hnl = 0.0_dp
755 IF (potential%sgp_pot%has_nonlocal) THEN
756 ALLOCATE (pgauss(1:m))
757 n = potential%sgp_pot%n_nonlocal
758 !
759 DO l = 0, potential%sgp_pot%lmax
760 cpassert(l <= ubound(basis%nbas, 1))
761 IF (.NOT. potential%sgp_pot%is_nonlocal(l)) cycle
762 ! overlap (a|p)
763 na = basis%nbas(l)
764 ALLOCATE (qmat(na, n))
765 DO i = 1, n
766 pgauss(:) = 0.0_dp
767 DO j = 1, n
768 a = potential%sgp_pot%a_nonlocal(j)
769 c = potential%sgp_pot%c_nonlocal(j, i, l)
770 pgauss(:) = pgauss(:) + c*exp(-a*rad(:)**2)*rad(:)**l
771 END DO
772 DO ia = 1, na
773 qmat(ia, i) = sum(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:))
774 END DO
775 END DO
776 DO i = 1, na
777 DO j = i, na
778 DO ia = 1, n
779 integrals%hnl(i, j, l) = integrals%hnl(i, j, l) &
780 + qmat(i, ia)*qmat(j, ia)*potential%sgp_pot%h_nonlocal(ia, l)
781 END DO
782 integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
783 END DO
784 END DO
785 DEALLOCATE (qmat)
786 END DO
787 DEALLOCATE (pgauss)
788 END IF
789
790 END SUBROUTINE sgpint_setup
791
792! **************************************************************************************************
793!> \brief ...
794!> \param integrals ...
795!> \param basis ...
796!> \param reltyp ...
797!> \param zcore ...
798!> \param alpha ...
799! **************************************************************************************************
800 SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
801 TYPE(atom_integrals), INTENT(INOUT) :: integrals
802 TYPE(atom_basis_type), INTENT(INOUT) :: basis
803 INTEGER, INTENT(IN) :: reltyp
804 REAL(dp), INTENT(IN) :: zcore
805 REAL(dp), INTENT(IN), OPTIONAL :: alpha
806
807 CHARACTER(len=*), PARAMETER :: routinen = 'atom_relint_setup'
808
809 INTEGER :: dkhorder, handle, i, k1, k2, l, m, n, nl
810 REAL(dp) :: ascal
811 REAL(dp), ALLOCATABLE, DIMENSION(:) :: cpot
812 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: modpot
813 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ener, sps
814 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hmat, pvp, sp, tp, vp, wfn
815
816 CALL timeset(routinen, handle)
817
818 SELECT CASE (reltyp)
819 CASE DEFAULT
820 cpabort("")
821 CASE (do_nonrel_atom, do_zoramp_atom, do_sczoramp_atom)
822 dkhorder = -1
823 CASE (do_dkh0_atom)
824 dkhorder = 0
825 CASE (do_dkh1_atom)
826 dkhorder = 1
827 CASE (do_dkh2_atom)
828 dkhorder = 2
829 CASE (do_dkh3_atom)
830 dkhorder = 3
831 END SELECT
832
833 SELECT CASE (reltyp)
834 CASE DEFAULT
835 cpabort("")
836 CASE (do_nonrel_atom)
837 ! nothing to do
838 NULLIFY (integrals%tzora, integrals%hdkh)
839 CASE (do_zoramp_atom, do_sczoramp_atom)
840
841 NULLIFY (integrals%hdkh)
842
843 IF (integrals%zorastat == 0) THEN
844 n = maxval(basis%nbas)
845 ALLOCATE (integrals%tzora(n, n, 0:lmat))
846 integrals%tzora = 0._dp
847 m = basis%grid%nr
848 ALLOCATE (modpot(1:m), cpot(1:m))
849 CALL calculate_model_potential(modpot, basis%grid, zcore)
850 ! Zora potential
851 cpot(1:m) = modpot(1:m)/(4._dp*c_light_au*c_light_au - 2._dp*modpot(1:m))
852 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
853 CALL numpot_matrix(integrals%tzora, cpot, basis, 0)
854 DO l = 0, lmat
855 nl = basis%nbas(l)
856 integrals%tzora(1:nl, 1:nl, l) = real(l*(l + 1), dp)*integrals%tzora(1:nl, 1:nl, l)
857 END DO
858 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
859 CALL numpot_matrix(integrals%tzora, cpot, basis, 2)
860 !
861 ! scaled ZORA
862 IF (reltyp == do_sczoramp_atom) THEN
863 ALLOCATE (hmat(n, n, 0:lmat), wfn(n, n, 0:lmat), ener(n, 0:lmat), pvp(n, n, 0:lmat), sps(n, n))
864 hmat(:, :, :) = integrals%kin + integrals%tzora
865 ! model potential
866 CALL numpot_matrix(hmat, modpot, basis, 0)
867 ! eigenvalues and eigenvectors
868 CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne, lmat)
869 ! relativistic kinetic energy
870 cpot(1:m) = c_light_au*c_light_au/(2._dp*c_light_au*c_light_au - modpot(1:m))**2
871 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
872 pvp = 0.0_dp
873 CALL numpot_matrix(pvp, cpot, basis, 0)
874 DO l = 0, lmat
875 nl = basis%nbas(l)
876 pvp(1:nl, 1:nl, l) = real(l*(l + 1), dp)*pvp(1:nl, 1:nl, l)
877 END DO
878 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
879 CALL numpot_matrix(pvp, cpot, basis, 2)
880 ! calculate psi*pvp*psi and the scaled orbital energies
881 ! actually, we directly calculate the energy difference
882 DO l = 0, lmat
883 nl = basis%nbas(l)
884 DO i = 1, integrals%nne(l)
885 IF (ener(i, l) < 0._dp) THEN
886 ascal = sum(wfn(1:nl, i, l)*matmul(pvp(1:nl, 1:nl, l), wfn(1:nl, i, l)))
887 ener(i, l) = ener(i, l)*ascal/(1.0_dp + ascal)
888 ELSE
889 ener(i, l) = 0.0_dp
890 END IF
891 END DO
892 END DO
893 ! correction term is calculated as a projector
894 hmat = 0.0_dp
895 DO l = 0, lmat
896 nl = basis%nbas(l)
897 DO i = 1, integrals%nne(l)
898 DO k1 = 1, nl
899 DO k2 = 1, nl
900 hmat(k1, k2, l) = hmat(k1, k2, l) + ener(i, l)*wfn(k1, i, l)*wfn(k2, i, l)
901 END DO
902 END DO
903 END DO
904 ! transform with overlap matrix
905 sps(1:nl, 1:nl) = matmul(integrals%ovlp(1:nl, 1:nl, l), &
906 matmul(hmat(1:nl, 1:nl, l), integrals%ovlp(1:nl, 1:nl, l)))
907 ! add scaling correction to tzora
908 integrals%tzora(1:nl, 1:nl, l) = integrals%tzora(1:nl, 1:nl, l) - sps(1:nl, 1:nl)
909 END DO
910
911 DEALLOCATE (hmat, wfn, ener, pvp, sps)
912 END IF
913 !
914 DEALLOCATE (modpot, cpot)
915
916 integrals%zorastat = 1
917
918 END IF
919
920 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
921
922 NULLIFY (integrals%tzora)
923
924 IF (integrals%dkhstat == 0) THEN
925 n = maxval(basis%nbas)
926 ALLOCATE (integrals%hdkh(n, n, 0:lmat))
927 integrals%hdkh = 0._dp
928
929 m = maxval(basis%nprim)
930 ALLOCATE (tp(m, m, 0:lmat), sp(m, m, 0:lmat), vp(m, m, 0:lmat), pvp(m, m, 0:lmat))
931 tp = 0._dp
932 sp = 0._dp
933 vp = 0._dp
934 pvp = 0._dp
935
936 SELECT CASE (basis%basis_type)
937 CASE DEFAULT
938 cpabort("")
939 CASE (gto_basis, cgto_basis)
940
941 DO l = 0, lmat
942 m = basis%nprim(l)
943 IF (m > 0) THEN
944 CALL sg_kinetic(tp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
945 CALL sg_overlap(sp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
946 IF (PRESENT(alpha)) THEN
947 CALL sg_erfc(vp(1:m, 1:m, l), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
948 ELSE
949 CALL sg_nuclear(vp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
950 END IF
951 CALL sg_kinnuc(pvp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
952 vp(1:m, 1:m, l) = -zcore*vp(1:m, 1:m, l)
953 pvp(1:m, 1:m, l) = -zcore*pvp(1:m, 1:m, l)
954 END IF
955 END DO
956
957 CASE (sto_basis)
958 cpabort("")
959 CASE (num_basis)
960 cpabort("")
961 END SELECT
962
963 CALL dkh_integrals(integrals, basis, dkhorder, sp, tp, vp, pvp)
964
965 integrals%dkhstat = 1
966
967 DEALLOCATE (tp, sp, vp, pvp)
968
969 ELSE
970 cpassert(ASSOCIATED(integrals%hdkh))
971 END IF
972
973 END SELECT
974
975 CALL timestop(handle)
976
977 END SUBROUTINE atom_relint_setup
978
979! **************************************************************************************************
980!> \brief ...
981!> \param integrals ...
982!> \param basis ...
983!> \param order ...
984!> \param sp ...
985!> \param tp ...
986!> \param vp ...
987!> \param pvp ...
988! **************************************************************************************************
989 SUBROUTINE dkh_integrals(integrals, basis, order, sp, tp, vp, pvp)
990 TYPE(atom_integrals), INTENT(INOUT) :: integrals
991 TYPE(atom_basis_type), INTENT(INOUT) :: basis
992 INTEGER, INTENT(IN) :: order
993 REAL(dp), DIMENSION(:, :, 0:) :: sp, tp, vp, pvp
994
995 INTEGER :: l, m, n
996 REAL(dp), DIMENSION(:, :, :), POINTER :: hdkh
997
998 cpassert(order >= 0)
999
1000 hdkh => integrals%hdkh
1001
1002 DO l = 0, lmat
1003 n = integrals%n(l)
1004 m = basis%nprim(l)
1005 IF (n > 0) THEN
1006 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)
1007 SELECT CASE (basis%basis_type)
1008 CASE DEFAULT
1009 cpabort("")
1010 CASE (gto_basis)
1011 cpassert(n == m)
1012 integrals%hdkh(1:n, 1:n, l) = tp(1:n, 1:n, l) + vp(1:n, 1:n, l)
1013 CASE (cgto_basis)
1014 CALL contract2(integrals%hdkh(1:n, 1:n, l), tp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1015 CALL contract2add(integrals%hdkh(1:n, 1:n, l), vp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1016 CASE (sto_basis)
1017 cpabort("")
1018 CASE (num_basis)
1019 cpabort("")
1020 END SELECT
1021 ELSE
1022 integrals%hdkh(1:n, 1:n, l) = 0._dp
1023 END IF
1024 END DO
1025
1026 END SUBROUTINE dkh_integrals
1027
1028! **************************************************************************************************
1029!> \brief Calculate overlap matrix between two atomic basis sets.
1030!> \param ovlap overlap matrix <chi_{a,l} | chi_{b,l}>
1031!> \param basis_a first basis set
1032!> \param basis_b second basis set
1033! **************************************************************************************************
1034 SUBROUTINE atom_basis_projection_overlap(ovlap, basis_a, basis_b)
1035 REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: ovlap
1036 TYPE(atom_basis_type), INTENT(IN) :: basis_a, basis_b
1037
1038 INTEGER :: i, j, l, ma, mb, na, nb
1039 LOGICAL :: ebas
1040 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
1041
1042 ebas = (basis_a%basis_type == basis_b%basis_type)
1043
1044 ovlap = 0.0_dp
1045
1046 IF (ebas) THEN
1047 SELECT CASE (basis_a%basis_type)
1048 CASE DEFAULT
1049 cpabort("")
1050 CASE (gto_basis)
1051 DO l = 0, lmat
1052 na = basis_a%nbas(l)
1053 nb = basis_b%nbas(l)
1054 CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l))
1055 END DO
1056 CASE (cgto_basis)
1057 DO l = 0, lmat
1058 na = basis_a%nbas(l)
1059 nb = basis_b%nbas(l)
1060 ma = basis_a%nprim(l)
1061 mb = basis_b%nprim(l)
1062 ALLOCATE (omat(ma, mb))
1063 CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l))
1064 ovlap(1:na, 1:nb, l) = matmul(transpose(basis_a%cm(1:ma, 1:na, l)), &
1065 matmul(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l)))
1066 DEALLOCATE (omat)
1067 END DO
1068 CASE (sto_basis)
1069 DO l = 0, lmat
1070 na = basis_a%nbas(l)
1071 nb = basis_b%nbas(l)
1072 CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), &
1073 basis_a%ns(1:na, l), basis_b%as(1:nb, l))
1074 END DO
1075 CASE (num_basis)
1076 cpassert(atom_compare_grids(basis_a%grid, basis_b%grid))
1077 DO l = 0, lmat
1078 na = basis_a%nbas(l)
1079 nb = basis_b%nbas(l)
1080 DO j = 1, nb
1081 DO i = 1, na
1082 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1083 END DO
1084 END DO
1085 END DO
1086 END SELECT
1087 ELSE
1088 cpassert(atom_compare_grids(basis_a%grid, basis_b%grid))
1089 DO l = 0, lmat
1090 na = basis_a%nbas(l)
1091 nb = basis_b%nbas(l)
1092 DO j = 1, nb
1093 DO i = 1, na
1094 ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1095 END DO
1096 END DO
1097 END DO
1098 END IF
1099
1100 END SUBROUTINE atom_basis_projection_overlap
1101
1102! **************************************************************************************************
1103!> \brief Release memory allocated for atomic integrals (valence electrons).
1104!> \param integrals atomic integrals
1105! **************************************************************************************************
1106 SUBROUTINE atom_int_release(integrals)
1107 TYPE(atom_integrals), INTENT(INOUT) :: integrals
1108
1109 INTEGER :: ll
1110
1111 IF (ASSOCIATED(integrals%ovlp)) THEN
1112 DEALLOCATE (integrals%ovlp)
1113 END IF
1114 IF (ASSOCIATED(integrals%kin)) THEN
1115 DEALLOCATE (integrals%kin)
1116 END IF
1117 IF (ASSOCIATED(integrals%conf)) THEN
1118 DEALLOCATE (integrals%conf)
1119 END IF
1120 DO ll = 1, SIZE(integrals%ceri)
1121 IF (ASSOCIATED(integrals%ceri(ll)%int)) THEN
1122 DEALLOCATE (integrals%ceri(ll)%int)
1123 END IF
1124 IF (ASSOCIATED(integrals%eeri(ll)%int)) THEN
1125 DEALLOCATE (integrals%eeri(ll)%int)
1126 END IF
1127 END DO
1128 IF (ASSOCIATED(integrals%utrans)) THEN
1129 DEALLOCATE (integrals%utrans)
1130 END IF
1131 IF (ASSOCIATED(integrals%uptrans)) THEN
1132 DEALLOCATE (integrals%uptrans)
1133 END IF
1134
1135 integrals%status = 0
1136
1137 END SUBROUTINE atom_int_release
1138
1139! **************************************************************************************************
1140!> \brief Release memory allocated for atomic integrals (core electrons).
1141!> \param integrals atomic integrals
1142! **************************************************************************************************
1143 SUBROUTINE atom_ppint_release(integrals)
1144 TYPE(atom_integrals), INTENT(INOUT) :: integrals
1145
1146 IF (ASSOCIATED(integrals%hnl)) THEN
1147 DEALLOCATE (integrals%hnl)
1148 END IF
1149 IF (ASSOCIATED(integrals%core)) THEN
1150 DEALLOCATE (integrals%core)
1151 END IF
1152 IF (ASSOCIATED(integrals%clsd)) THEN
1153 DEALLOCATE (integrals%clsd)
1154 END IF
1155
1156 integrals%ppstat = 0
1157
1158 END SUBROUTINE atom_ppint_release
1159
1160! **************************************************************************************************
1161!> \brief Release memory allocated for atomic integrals (relativistic effects).
1162!> \param integrals atomic integrals
1163! **************************************************************************************************
1164 SUBROUTINE atom_relint_release(integrals)
1165 TYPE(atom_integrals), INTENT(INOUT) :: integrals
1166
1167 IF (ASSOCIATED(integrals%tzora)) THEN
1168 DEALLOCATE (integrals%tzora)
1169 END IF
1170 IF (ASSOCIATED(integrals%hdkh)) THEN
1171 DEALLOCATE (integrals%hdkh)
1172 END IF
1173
1174 integrals%zorastat = 0
1175 integrals%dkhstat = 0
1176
1177 END SUBROUTINE atom_relint_release
1178
1179! **************************************************************************************************
1180!> \brief Calculate model potential. It is useful to guess atomic orbitals.
1181!> \param modpot model potential
1182!> \param grid atomic radial grid
1183!> \param zcore nuclear charge
1184! **************************************************************************************************
1185 SUBROUTINE calculate_model_potential(modpot, grid, zcore)
1186 REAL(dp), DIMENSION(:), INTENT(INOUT) :: modpot
1187 TYPE(grid_atom_type), INTENT(IN) :: grid
1188 REAL(dp), INTENT(IN) :: zcore
1189
1190 INTEGER :: i, ii, l, ll, n, z
1191 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: pot, rho
1192 TYPE(atom_state) :: state
1193
1194 n = SIZE(modpot)
1195 ALLOCATE (rho(n), pot(n))
1196
1197 ! fill default occupation
1198 state%core = 0._dp
1199 state%occ = 0._dp
1200 state%multiplicity = -1
1201 z = nint(zcore)
1202 DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
1203 IF (ptable(z)%e_conv(l) /= 0) THEN
1204 state%maxl_occ = l
1205 ll = 2*(2*l + 1)
1206 DO i = 1, SIZE(state%occ, 2)
1207 ii = ptable(z)%e_conv(l) - (i - 1)*ll
1208 IF (ii <= ll) THEN
1209 state%occ(l, i) = ii
1210 EXIT
1211 ELSE
1212 state%occ(l, i) = ll
1213 END IF
1214 END DO
1215 END IF
1216 END DO
1217
1218 modpot = -zcore/grid%rad(:)
1219
1220 ! Coulomb potential
1221 CALL slater_density(rho, pot, nint(zcore), state, grid)
1222 CALL coulomb_potential_numeric(pot, rho, grid)
1223 modpot = modpot + pot
1224
1225 ! XC potential
1226 CALL wigner_slater_functional(rho, pot)
1227 modpot = modpot + pot
1228
1229 DEALLOCATE (rho, pot)
1230
1231 END SUBROUTINE calculate_model_potential
1232
1233END 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:1144
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:801
atom_operators::atom_relint_release
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
Definition atom_operators.F:1165
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:1035
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:1186
atom_operators::atom_int_release
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Definition atom_operators.F:1107
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:3023
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: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: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