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