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