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 11154 : 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 11154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, w, work
75 11154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat, vmat
76 11154 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: eri
77 :
78 11154 : CALL timeset(routineN, handle)
79 :
80 11154 : IF (integrals%status == 0) THEN
81 78050 : n = MAXVAL(basis%nbas)
82 78050 : integrals%n = basis%nbas
83 :
84 11150 : IF (PRESENT(eri_coulomb)) THEN
85 11122 : integrals%eri_coulomb = eri_coulomb
86 : ELSE
87 28 : integrals%eri_coulomb = .FALSE.
88 : END IF
89 11150 : IF (PRESENT(eri_exchange)) THEN
90 11124 : integrals%eri_exchange = eri_exchange
91 : ELSE
92 26 : integrals%eri_exchange = .FALSE.
93 : END IF
94 11150 : IF (PRESENT(all_nu)) THEN
95 0 : integrals%all_nu = all_nu
96 : ELSE
97 11150 : integrals%all_nu = .FALSE.
98 : END IF
99 :
100 11150 : NULLIFY (integrals%ovlp, integrals%kin, integrals%core, integrals%conf)
101 1126150 : DO ll = 1, SIZE(integrals%ceri)
102 1126150 : NULLIFY (integrals%ceri(ll)%int, integrals%eeri(ll)%int)
103 : END DO
104 :
105 55738 : ALLOCATE (integrals%ovlp(n, n, 0:lmat))
106 2683574 : integrals%ovlp = 0._dp
107 :
108 33438 : ALLOCATE (integrals%kin(n, n, 0:lmat))
109 2683574 : integrals%kin = 0._dp
110 :
111 11150 : integrals%status = 1
112 :
113 11150 : IF (PRESENT(potential)) THEN
114 11122 : IF (potential%confinement) THEN
115 25744 : ALLOCATE (integrals%conf(n, n, 0:lmat))
116 1060456 : integrals%conf = 0._dp
117 8584 : m = basis%grid%nr
118 25752 : ALLOCATE (cpot(1:m))
119 8584 : IF (potential%conf_type == poly_conf) THEN
120 8364 : rc = potential%rcon
121 8364 : sc = potential%scon
122 3351132 : 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 8584 : CALL numpot_matrix(integrals%conf, cpot, basis, 0)
144 8584 : DEALLOCATE (cpot)
145 : END IF
146 : END IF
147 :
148 11150 : SELECT CASE (basis%basis_type)
149 : CASE DEFAULT
150 0 : CPABORT("")
151 : CASE (GTO_BASIS)
152 9800 : DO l = 0, lmat
153 8400 : n = integrals%n(l)
154 8400 : CALL sg_overlap(integrals%ovlp(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
155 9800 : 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 1400 : 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 1400 : 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 60340 : DO l = 0, lmat
202 51720 : n = integrals%n(l)
203 51720 : m = basis%nprim(l)
204 60340 : IF (n > 0 .AND. m > 0) THEN
205 79204 : ALLOCATE (omat(m, m))
206 :
207 19801 : CALL sg_overlap(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
208 19801 : CALL contract2(integrals%ovlp(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
209 19801 : CALL sg_kinetic(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
210 19801 : CALL contract2(integrals%kin(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
211 19801 : DEALLOCATE (omat)
212 : END IF
213 : END DO
214 8620 : 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 8620 : 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 11150 : CPABORT("")
287 : END SELECT
288 :
289 : ! setup transformation matrix to get an orthogonal basis, remove linear dependencies
290 11150 : NULLIFY (integrals%utrans, integrals%uptrans)
291 78050 : n = MAXVAL(basis%nbas)
292 78026 : ALLOCATE (integrals%utrans(n, n, 0:lmat), integrals%uptrans(n, n, 0:lmat))
293 2683574 : integrals%utrans = 0._dp
294 2683574 : integrals%uptrans = 0._dp
295 78050 : integrals%nne = integrals%n
296 11150 : lwork = 10*n
297 111464 : ALLOCATE (omat(n, n), vmat(n, n), w(n), work(lwork))
298 78050 : DO l = 0, lmat
299 66900 : n = integrals%n(l)
300 78050 : IF (n > 0) THEN
301 1716521 : omat(1:n, 1:n) = integrals%ovlp(1:n, 1:n, l)
302 25809 : CALL jacobi(omat(1:n, 1:n), w(1:n), vmat(1:n, 1:n))
303 1716521 : omat(1:n, 1:n) = vmat(1:n, 1:n)
304 25809 : ii = 0
305 124288 : DO i = 1, n
306 124288 : IF (w(i) > basis%eps_eig) THEN
307 83449 : ii = ii + 1
308 1102592 : integrals%utrans(1:n, ii, l) = omat(1:n, i)/SQRT(w(i))
309 : END IF
310 : END DO
311 25809 : integrals%nne(l) = ii
312 25809 : IF (ii > 0) THEN
313 15095552 : omat(1:ii, 1:ii) = MATMUL(TRANSPOSE(integrals%utrans(1:n, 1:ii, l)), integrals%utrans(1:n, 1:ii, l))
314 109258 : DO i = 1, ii
315 109258 : integrals%uptrans(i, i, l) = 1._dp
316 : END DO
317 1988641 : CALL dgesv(ii, ii, omat(1:ii, 1:ii), ii, ipiv, integrals%uptrans(1:ii, 1:ii, l), ii, info)
318 25809 : CPASSERT(info == 0)
319 : END IF
320 : END IF
321 : END DO
322 11150 : DEALLOCATE (omat, vmat, w, work)
323 : END IF
324 :
325 11154 : CALL timestop(handle)
326 :
327 22308 : END SUBROUTINE atom_int_setup
328 :
329 : ! **************************************************************************************************
330 : !> \brief ...
331 : !> \param integrals ...
332 : !> \param basis ...
333 : !> \param potential ...
334 : ! **************************************************************************************************
335 11484 : 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 11484 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, xmat
345 11484 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat, spmat
346 11484 : REAL(KIND=dp), DIMENSION(:), POINTER :: rad
347 :
348 11484 : CALL timeset(routineN, handle)
349 :
350 11484 : IF (integrals%ppstat == 0) THEN
351 80360 : n = MAXVAL(basis%nbas)
352 80360 : integrals%n = basis%nbas
353 :
354 11480 : NULLIFY (integrals%core, integrals%hnl)
355 :
356 57388 : ALLOCATE (integrals%hnl(n, n, 0:lmat))
357 3441080 : integrals%hnl = 0._dp
358 :
359 34428 : ALLOCATE (integrals%core(n, n, 0:lmat))
360 3441080 : integrals%core = 0._dp
361 :
362 34428 : ALLOCATE (integrals%clsd(n, n, 0:lmat))
363 3441080 : integrals%clsd = 0._dp
364 :
365 11480 : integrals%ppstat = 1
366 :
367 11480 : SELECT CASE (basis%basis_type)
368 : CASE DEFAULT
369 0 : CPABORT("")
370 : CASE (GTO_BASIS)
371 :
372 10426 : 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 1336 : alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp)
380 9352 : DO l = 0, lmat
381 8016 : n = integrals%n(l)
382 37272 : ALLOCATE (omat(n, n), spmat(n, 5))
383 :
384 1045340 : omat = 0._dp
385 8016 : CALL sg_erf(omat(1:n, 1:n), l, alpha, basis%am(1:n, l), basis%am(1:n, l))
386 1045340 : integrals%core(1:n, 1:n, l) = -potential%gth_pot%zion*omat(1:n, 1:n)
387 18720 : DO i = 1, potential%gth_pot%ncl
388 1982092 : omat = 0._dp
389 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))
390 : integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
391 1990108 : potential%gth_pot%cl(i)*omat(1:n, 1:n)
392 : END DO
393 8016 : IF (potential%gth_pot%lpotextended) THEN
394 96 : DO k = 1, potential%gth_pot%nexp_lpot
395 228 : DO i = 1, potential%gth_pot%nct_lpot(k)
396 36036 : omat = 0._dp
397 : CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lpot(k), l, &
398 132 : basis%am(1:n, l), basis%am(1:n, l))
399 : integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
400 36096 : potential%gth_pot%cval_lpot(i, k)*omat(1:n, 1:n)
401 : END DO
402 : END DO
403 : END IF
404 8016 : IF (potential%gth_pot%lsdpot) THEN
405 0 : DO k = 1, potential%gth_pot%nexp_lsd
406 0 : DO i = 1, potential%gth_pot%nct_lsd(k)
407 0 : omat = 0._dp
408 : CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lsd(k), l, &
409 0 : basis%am(1:n, l), basis%am(1:n, l))
410 : integrals%clsd(1:n, 1:n, l) = integrals%clsd(1:n, 1:n, l) + &
411 0 : potential%gth_pot%cval_lsd(i, k)*omat(1:n, 1:n)
412 : END DO
413 : END DO
414 : END IF
415 :
416 303096 : spmat = 0._dp
417 8016 : m = potential%gth_pot%nl(l)
418 13844 : DO i = 1, m
419 13844 : CALL sg_proj_ol(spmat(1:n, i), l, basis%am(1:n, l), i - 1, potential%gth_pot%rcnl(l))
420 : END DO
421 8016 : integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:m), &
422 2355606 : MATMUL(potential%gth_pot%hnl(1:m, 1:m, l), TRANSPOSE(spmat(1:n, 1:m))))
423 :
424 9352 : DEALLOCATE (omat, spmat)
425 : END DO
426 : CASE (upf_pseudo)
427 4 : CALL upfint_setup(integrals, basis, potential)
428 : CASE (sgp_pseudo)
429 0 : CALL sgpint_setup(integrals, basis, potential)
430 : CASE DEFAULT
431 1574 : CPABORT("")
432 : END SELECT
433 :
434 : CASE (CGTO_BASIS)
435 :
436 11090 : SELECT CASE (potential%ppot_type)
437 : CASE (no_pseudo, ecp_pseudo)
438 8288 : DO l = 0, lmat
439 7104 : n = integrals%n(l)
440 7104 : m = basis%nprim(l)
441 19596 : ALLOCATE (omat(m, m))
442 :
443 7104 : CALL sg_nuclear(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
444 7104 : CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
445 :
446 8288 : DEALLOCATE (omat)
447 : END DO
448 : CASE (gth_pseudo)
449 7422 : alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp)
450 51954 : DO l = 0, lmat
451 44532 : n = integrals%n(l)
452 44532 : m = basis%nprim(l)
453 51954 : IF (n > 0 .AND. m > 0) THEN
454 136568 : ALLOCATE (omat(m, m), spmat(n, 5), xmat(m))
455 :
456 376615 : omat = 0._dp
457 17071 : CALL sg_erf(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
458 376615 : omat(1:m, 1:m) = -potential%gth_pot%zion*omat(1:m, 1:m)
459 17071 : CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
460 49995 : DO i = 1, potential%gth_pot%ncl
461 722396 : omat = 0._dp
462 32924 : CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%rc, l, basis%am(1:m, l), basis%am(1:m, l))
463 722396 : omat(1:m, 1:m) = potential%gth_pot%cl(i)*omat(1:m, 1:m)
464 49995 : CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
465 : END DO
466 17071 : IF (potential%gth_pot%lpotextended) THEN
467 72 : DO k = 1, potential%gth_pot%nexp_lpot
468 168 : DO i = 1, potential%gth_pot%nct_lpot(k)
469 3128 : omat = 0._dp
470 : CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lpot(k), l, &
471 96 : basis%am(1:m, l), basis%am(1:m, l))
472 3128 : omat(1:m, 1:m) = potential%gth_pot%cval_lpot(i, k)*omat(1:m, 1:m)
473 140 : CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
474 : END DO
475 : END DO
476 : END IF
477 17071 : IF (potential%gth_pot%lsdpot) THEN
478 0 : DO k = 1, potential%gth_pot%nexp_lsd
479 0 : DO i = 1, potential%gth_pot%nct_lsd(k)
480 0 : omat = 0._dp
481 : CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lsd(k), l, &
482 0 : basis%am(1:m, l), basis%am(1:m, l))
483 0 : omat(1:m, 1:m) = potential%gth_pot%cval_lsd(i, k)*omat(1:m, 1:m)
484 0 : CALL contract2add(integrals%clsd(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
485 : END DO
486 : END DO
487 : END IF
488 :
489 243866 : spmat = 0._dp
490 17071 : k = potential%gth_pot%nl(l)
491 22376 : DO i = 1, k
492 5305 : CALL sg_proj_ol(xmat(1:m), l, basis%am(1:m, l), i - 1, potential%gth_pot%rcnl(l))
493 22376 : spmat(1:n, i) = MATMUL(TRANSPOSE(basis%cm(1:m, 1:n, l)), xmat(1:m))
494 : END DO
495 17071 : IF (k > 0) THEN
496 4563 : integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:k), &
497 18252 : MATMUL(potential%gth_pot%hnl(1:k, 1:k, l), &
498 161125 : TRANSPOSE(spmat(1:n, 1:k))))
499 : END IF
500 17071 : DEALLOCATE (omat, spmat, xmat)
501 : END IF
502 : END DO
503 : CASE (upf_pseudo)
504 0 : CALL upfint_setup(integrals, basis, potential)
505 : CASE (sgp_pseudo)
506 12 : CALL sgpint_setup(integrals, basis, potential)
507 : CASE DEFAULT
508 8618 : CPABORT("")
509 : END SELECT
510 :
511 : CASE (STO_BASIS)
512 :
513 2336 : SELECT CASE (potential%ppot_type)
514 : CASE (no_pseudo, ecp_pseudo)
515 7336 : DO l = 0, lmat
516 6288 : n = integrals%n(l)
517 : CALL sto_nuclear(integrals%core(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
518 7336 : basis%ns(1:n, l), basis%as(1:n, l))
519 : END DO
520 : CASE (gth_pseudo)
521 240 : rad => basis%grid%rad
522 240 : m = basis%grid%nr
523 720 : ALLOCATE (cpot(1:m))
524 240 : rc = potential%gth_pot%rc
525 240 : alpha = 1._dp/rc/SQRT(2._dp)
526 :
527 : ! local pseudopotential, we use erf = 1/r - erfc
528 5040 : integrals%core = 0._dp
529 102640 : DO i = 1, m
530 102640 : cpot(i) = potential%gth_pot%zion*erfc(alpha*rad(i))/rad(i)
531 : END DO
532 720 : DO i = 1, potential%gth_pot%ncl
533 480 : ii = 2*(i - 1)
534 205520 : cpot(1:m) = cpot(1:m) + potential%gth_pot%cl(i)*(rad/rc)**ii*EXP(-0.5_dp*(rad/rc)**2)
535 : END DO
536 240 : IF (potential%gth_pot%lpotextended) THEN
537 0 : DO k = 1, potential%gth_pot%nexp_lpot
538 0 : al = potential%gth_pot%alpha_lpot(k)
539 0 : DO i = 1, potential%gth_pot%nct_lpot(k)
540 0 : ii = 2*(i - 1)
541 0 : cpot(1:m) = cpot(1:m) + potential%gth_pot%cval_lpot(i, k)*(rad/al)**ii*EXP(-0.5_dp*(rad/al)**2)
542 : END DO
543 : END DO
544 : END IF
545 240 : CALL numpot_matrix(integrals%core, cpot, basis, 0)
546 1680 : DO l = 0, lmat
547 1440 : n = integrals%n(l)
548 3560 : ALLOCATE (omat(n, n))
549 2280 : omat = 0._dp
550 : CALL sto_nuclear(omat(1:n, 1:n), basis%ns(1:n, l), basis%as(1:n, l), &
551 1440 : basis%ns(1:n, l), basis%as(1:n, l))
552 2280 : integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) - potential%gth_pot%zion*omat(1:n, 1:n)
553 1680 : DEALLOCATE (omat)
554 : END DO
555 :
556 240 : IF (potential%gth_pot%lsdpot) THEN
557 0 : cpot = 0._dp
558 0 : DO k = 1, potential%gth_pot%nexp_lsd
559 0 : al = potential%gth_pot%alpha_lsd(k)
560 0 : DO i = 1, potential%gth_pot%nct_lsd(k)
561 0 : ii = 2*(i - 1)
562 0 : cpot(:) = cpot + potential%gth_pot%cval_lsd(i, k)*(rad/al)**ii*EXP(-0.5_dp*(rad/al)**2)
563 : END DO
564 : END DO
565 0 : CALL numpot_matrix(integrals%clsd, cpot, basis, 0)
566 : END IF
567 :
568 1680 : DO l = 0, lmat
569 1440 : n = integrals%n(l)
570 : ! non local pseudopotential
571 3220 : ALLOCATE (spmat(n, 5))
572 10440 : spmat = 0._dp
573 1440 : k = potential%gth_pot%nl(l)
574 1688 : DO i = 1, k
575 248 : rc = potential%gth_pot%rcnl(l)
576 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)
577 1872 : DO j = 1, basis%nbas(l)
578 432 : spmat(j, i) = integrate_grid(cpot, basis%bf(:, j, l), basis%grid)
579 : END DO
580 : END DO
581 1440 : integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:k), &
582 3468 : MATMUL(potential%gth_pot%hnl(1:k, 1:k, l), &
583 4768 : TRANSPOSE(spmat(1:n, 1:k))))
584 1680 : DEALLOCATE (spmat)
585 : END DO
586 240 : DEALLOCATE (cpot)
587 :
588 : CASE (upf_pseudo)
589 0 : CALL upfint_setup(integrals, basis, potential)
590 : CASE (sgp_pseudo)
591 0 : CALL sgpint_setup(integrals, basis, potential)
592 : CASE DEFAULT
593 1288 : CPABORT("")
594 : END SELECT
595 :
596 : CASE (NUM_BASIS)
597 11480 : CPABORT("")
598 : END SELECT
599 :
600 : ! add ecp_pseudo using numerical representation of basis
601 11480 : IF (potential%ppot_type == ecp_pseudo) THEN
602 : ! scale 1/r potential
603 10858 : integrals%core = -potential%ecp_pot%zion*integrals%core
604 : ! local potential
605 22 : m = basis%grid%nr
606 22 : rad => basis%grid%rad
607 66 : ALLOCATE (cpot(1:m))
608 8822 : cpot = 0._dp
609 76 : DO k = 1, potential%ecp_pot%nloc
610 54 : n = potential%ecp_pot%nrloc(k)
611 54 : alpha = potential%ecp_pot%bloc(k)
612 21676 : cpot(:) = cpot + potential%ecp_pot%aloc(k)*rad**(n - 2)*EXP(-alpha*rad**2)
613 : END DO
614 22 : CALL numpot_matrix(integrals%core, cpot, basis, 0)
615 : ! non local pseudopotential
616 44 : DO l = 0, MIN(potential%ecp_pot%lmax, lmat)
617 8822 : cpot = 0._dp
618 52 : DO k = 1, potential%ecp_pot%npot(l)
619 30 : n = potential%ecp_pot%nrpot(k, l)
620 30 : alpha = potential%ecp_pot%bpot(k, l)
621 12052 : cpot(:) = cpot + potential%ecp_pot%apot(k, l)*rad**(n - 2)*EXP(-alpha*rad**2)
622 : END DO
623 396 : DO i = 1, basis%nbas(l)
624 3366 : DO j = i, basis%nbas(l)
625 2992 : integrals%hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
626 3344 : integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
627 : END DO
628 : END DO
629 : END DO
630 22 : DEALLOCATE (cpot)
631 : END IF
632 :
633 : END IF
634 :
635 11484 : CALL timestop(handle)
636 :
637 22968 : END SUBROUTINE atom_ppint_setup
638 :
639 : ! **************************************************************************************************
640 : !> \brief ...
641 : !> \param integrals ...
642 : !> \param basis ...
643 : !> \param potential ...
644 : ! **************************************************************************************************
645 4 : SUBROUTINE upfint_setup(integrals, basis, potential)
646 : TYPE(atom_integrals), INTENT(INOUT) :: integrals
647 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
648 : TYPE(atom_potential_type), INTENT(IN) :: potential
649 :
650 : CHARACTER(len=4) :: ptype
651 : INTEGER :: i, j, k1, k2, la, lb, m, n
652 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: spot
653 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: spmat
654 : TYPE(atom_basis_type) :: gbasis
655 :
656 : ! get basis representation on UPF grid
657 4 : CALL atom_basis_gridrep(basis, gbasis, potential%upf_pot%r, potential%upf_pot%rab)
658 :
659 : ! local pseudopotential
660 6556 : integrals%core = 0._dp
661 4 : CALL numpot_matrix(integrals%core, potential%upf_pot%vlocal, gbasis, 0)
662 :
663 4 : ptype = ADJUSTL(TRIM(potential%upf_pot%pseudo_type))
664 4 : IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
665 : ! non local pseudopotential
666 14 : n = MAXVAL(integrals%n(:))
667 2 : m = potential%upf_pot%number_of_proj
668 8 : ALLOCATE (spmat(n, m))
669 36 : spmat = 0.0_dp
670 4 : DO i = 1, m
671 2 : la = potential%upf_pot%lbeta(i)
672 36 : DO j = 1, gbasis%nbas(la)
673 34 : spmat(j, i) = integrate_grid(potential%upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
674 : END DO
675 : END DO
676 4 : DO i = 1, m
677 2 : la = potential%upf_pot%lbeta(i)
678 6 : DO j = 1, m
679 2 : lb = potential%upf_pot%lbeta(j)
680 4 : IF (la == lb) THEN
681 34 : DO k1 = 1, gbasis%nbas(la)
682 546 : DO k2 = 1, gbasis%nbas(la)
683 : integrals%hnl(k1, k2, la) = integrals%hnl(k1, k2, la) + &
684 544 : spmat(k1, i)*potential%upf_pot%dion(i, j)*spmat(k2, j)
685 : END DO
686 : END DO
687 : END IF
688 : END DO
689 : END DO
690 2 : DEALLOCATE (spmat)
691 2 : ELSE IF (ptype(1:2) == "SL") THEN
692 : ! semi local pseudopotential
693 10 : DO la = 0, potential%upf_pot%l_max
694 8 : IF (la == potential%upf_pot%l_local) CYCLE
695 6 : m = SIZE(potential%upf_pot%vsemi(:, la + 1))
696 18 : ALLOCATE (spot(m))
697 2772 : spot(:) = potential%upf_pot%vsemi(:, la + 1) - potential%upf_pot%vlocal(:)
698 6 : n = basis%nbas(la)
699 102 : DO i = 1, n
700 918 : DO j = i, n
701 : integrals%core(i, j, la) = integrals%core(i, j, la) + &
702 : integrate_grid(spot(:), &
703 816 : gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
704 912 : integrals%core(j, i, la) = integrals%core(i, j, la)
705 : END DO
706 : END DO
707 10 : DEALLOCATE (spot)
708 : END DO
709 : ELSE
710 0 : CPABORT("Pseudopotential type: ["//ADJUSTL(TRIM(ptype))//"] not known")
711 : END IF
712 :
713 : ! release basis representation on UPF grid
714 4 : CALL release_atom_basis(gbasis)
715 :
716 80 : END SUBROUTINE upfint_setup
717 :
718 : ! **************************************************************************************************
719 : !> \brief ...
720 : !> \param integrals ...
721 : !> \param basis ...
722 : !> \param potential ...
723 : ! **************************************************************************************************
724 12 : SUBROUTINE sgpint_setup(integrals, basis, potential)
725 : TYPE(atom_integrals), INTENT(INOUT) :: integrals
726 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
727 : TYPE(atom_potential_type), INTENT(IN) :: potential
728 :
729 : INTEGER :: i, ia, j, l, m, n, na
730 : REAL(KIND=dp) :: a, c, rc, zval
731 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, pgauss
732 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qmat
733 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: rad
734 :
735 12 : rad => basis%grid%rad
736 12 : m = basis%grid%nr
737 :
738 : ! local pseudopotential
739 1524 : integrals%core = 0._dp
740 36 : ALLOCATE (cpot(m))
741 4812 : cpot = 0.0_dp
742 12 : zval = potential%sgp_pot%zion
743 4812 : DO i = 1, m
744 4800 : rc = rad(i)/potential%sgp_pot%ac_local/SQRT(2.0_dp)
745 4812 : cpot(i) = cpot(i) - zval/rad(i)*erf(rc)
746 : END DO
747 156 : DO i = 1, potential%sgp_pot%n_local
748 57756 : cpot(:) = cpot(:) + potential%sgp_pot%c_local(i)*EXP(-potential%sgp_pot%a_local(i)*rad(:)**2)
749 : END DO
750 12 : CALL numpot_matrix(integrals%core, cpot, basis, 0)
751 12 : DEALLOCATE (cpot)
752 :
753 : ! nonlocal pseudopotential
754 1524 : integrals%hnl = 0.0_dp
755 12 : IF (potential%sgp_pot%has_nonlocal) THEN
756 18 : ALLOCATE (pgauss(1:m))
757 6 : n = potential%sgp_pot%n_nonlocal
758 : !
759 12 : DO l = 0, potential%sgp_pot%lmax
760 6 : CPASSERT(l <= UBOUND(basis%nbas, 1))
761 6 : IF (.NOT. potential%sgp_pot%is_nonlocal(l)) CYCLE
762 : ! overlap (a|p)
763 6 : na = basis%nbas(l)
764 24 : ALLOCATE (qmat(na, n))
765 54 : DO i = 1, n
766 19248 : pgauss(:) = 0.0_dp
767 432 : DO j = 1, n
768 384 : a = potential%sgp_pot%a_nonlocal(j)
769 384 : c = potential%sgp_pot%c_nonlocal(j, i, l)
770 154032 : pgauss(:) = pgauss(:) + c*EXP(-a*rad(:)**2)*rad(:)**l
771 : END DO
772 246 : DO ia = 1, na
773 77040 : qmat(ia, i) = SUM(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:))
774 : END DO
775 : END DO
776 30 : DO i = 1, na
777 90 : DO j = i, na
778 540 : DO ia = 1, n
779 : integrals%hnl(i, j, l) = integrals%hnl(i, j, l) &
780 540 : + qmat(i, ia)*qmat(j, ia)*potential%sgp_pot%h_nonlocal(ia, l)
781 : END DO
782 84 : integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
783 : END DO
784 : END DO
785 12 : DEALLOCATE (qmat)
786 : END DO
787 6 : DEALLOCATE (pgauss)
788 : END IF
789 :
790 24 : END SUBROUTINE sgpint_setup
791 :
792 : ! **************************************************************************************************
793 : !> \brief ...
794 : !> \param integrals ...
795 : !> \param basis ...
796 : !> \param reltyp ...
797 : !> \param zcore ...
798 : !> \param alpha ...
799 : ! **************************************************************************************************
800 9908 : SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
801 : TYPE(atom_integrals), INTENT(INOUT) :: integrals
802 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
803 : INTEGER, INTENT(IN) :: reltyp
804 : REAL(dp), INTENT(IN) :: zcore
805 : REAL(dp), INTENT(IN), OPTIONAL :: alpha
806 :
807 : CHARACTER(len=*), PARAMETER :: routineN = 'atom_relint_setup'
808 :
809 : INTEGER :: dkhorder, handle, i, k1, k2, l, m, n, nl
810 : REAL(dp) :: ascal
811 9908 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: cpot
812 9908 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: modpot
813 9908 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ener, sps
814 9908 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hmat, pvp, sp, tp, vp, wfn
815 :
816 9908 : CALL timeset(routineN, handle)
817 :
818 9908 : SELECT CASE (reltyp)
819 : CASE DEFAULT
820 0 : CPABORT("")
821 : CASE (do_nonrel_atom, do_zoramp_atom, do_sczoramp_atom)
822 9836 : dkhorder = -1
823 : CASE (do_dkh0_atom)
824 2 : dkhorder = 0
825 : CASE (do_dkh1_atom)
826 2 : dkhorder = 1
827 : CASE (do_dkh2_atom)
828 30 : dkhorder = 2
829 : CASE (do_dkh3_atom)
830 9908 : dkhorder = 3
831 : END SELECT
832 :
833 0 : SELECT CASE (reltyp)
834 : CASE DEFAULT
835 0 : CPABORT("")
836 : CASE (do_nonrel_atom)
837 : ! nothing to do
838 9808 : NULLIFY (integrals%tzora, integrals%hdkh)
839 : CASE (do_zoramp_atom, do_sczoramp_atom)
840 :
841 28 : NULLIFY (integrals%hdkh)
842 :
843 28 : IF (integrals%zorastat == 0) THEN
844 196 : n = MAXVAL(basis%nbas)
845 140 : ALLOCATE (integrals%tzora(n, n, 0:lmat))
846 133636 : integrals%tzora = 0._dp
847 28 : m = basis%grid%nr
848 112 : ALLOCATE (modpot(1:m), cpot(1:m))
849 28 : CALL calculate_model_potential(modpot, basis%grid, zcore)
850 : ! Zora potential
851 11228 : cpot(1:m) = modpot(1:m)/(4._dp*c_light_au*c_light_au - 2._dp*modpot(1:m))
852 11228 : cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
853 28 : CALL numpot_matrix(integrals%tzora, cpot, basis, 0)
854 196 : DO l = 0, lmat
855 168 : nl = basis%nbas(l)
856 123660 : integrals%tzora(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*integrals%tzora(1:nl, 1:nl, l)
857 : END DO
858 11228 : cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
859 28 : CALL numpot_matrix(integrals%tzora, cpot, basis, 2)
860 : !
861 : ! scaled ZORA
862 28 : IF (reltyp == do_sczoramp_atom) THEN
863 168 : ALLOCATE (hmat(n, n, 0:lmat), wfn(n, n, 0:lmat), ener(n, 0:lmat), pvp(n, n, 0:lmat), sps(n, n))
864 31946 : hmat(:, :, :) = integrals%kin + integrals%tzora
865 : ! model potential
866 14 : CALL numpot_matrix(hmat, modpot, basis, 0)
867 : ! eigenvalues and eigenvectors
868 14 : CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne, lmat)
869 : ! relativistic kinetic energy
870 5614 : cpot(1:m) = c_light_au*c_light_au/(2._dp*c_light_au*c_light_au - modpot(1:m))**2
871 5614 : cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
872 31946 : pvp = 0.0_dp
873 14 : CALL numpot_matrix(pvp, cpot, basis, 0)
874 98 : DO l = 0, lmat
875 84 : nl = basis%nbas(l)
876 24010 : pvp(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*pvp(1:nl, 1:nl, l)
877 : END DO
878 5614 : cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
879 14 : CALL numpot_matrix(pvp, cpot, basis, 2)
880 : ! calculate psi*pvp*psi and the scaled orbital energies
881 : ! actually, we directly calculate the energy difference
882 98 : DO l = 0, lmat
883 84 : nl = basis%nbas(l)
884 578 : DO i = 1, integrals%nne(l)
885 564 : IF (ener(i, l) < 0._dp) THEN
886 21276 : ascal = SUM(wfn(1:nl, i, l)*MATMUL(pvp(1:nl, 1:nl, l), wfn(1:nl, i, l)))
887 28 : ener(i, l) = ener(i, l)*ascal/(1.0_dp + ascal)
888 : ELSE
889 452 : ener(i, l) = 0.0_dp
890 : END IF
891 : END DO
892 : END DO
893 : ! correction term is calculated as a projector
894 31946 : hmat = 0.0_dp
895 98 : DO l = 0, lmat
896 84 : nl = basis%nbas(l)
897 564 : DO i = 1, integrals%nne(l)
898 14238 : DO k1 = 1, nl
899 494628 : DO k2 = 1, nl
900 494148 : hmat(k1, k2, l) = hmat(k1, k2, l) + ener(i, l)*wfn(k1, i, l)*wfn(k2, i, l)
901 : END DO
902 : END DO
903 : END DO
904 : ! transform with overlap matrix
905 84 : sps(1:nl, 1:nl) = MATMUL(integrals%ovlp(1:nl, 1:nl, l), &
906 410216 : MATMUL(hmat(1:nl, 1:nl, l), integrals%ovlp(1:nl, 1:nl, l)))
907 : ! add scaling correction to tzora
908 24010 : integrals%tzora(1:nl, 1:nl, l) = integrals%tzora(1:nl, 1:nl, l) - sps(1:nl, 1:nl)
909 : END DO
910 :
911 14 : DEALLOCATE (hmat, wfn, ener, pvp, sps)
912 : END IF
913 : !
914 28 : DEALLOCATE (modpot, cpot)
915 :
916 28 : integrals%zorastat = 1
917 :
918 : END IF
919 :
920 : CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
921 :
922 72 : NULLIFY (integrals%tzora)
923 :
924 9908 : IF (integrals%dkhstat == 0) THEN
925 504 : n = MAXVAL(basis%nbas)
926 360 : ALLOCATE (integrals%hdkh(n, n, 0:lmat))
927 321120 : integrals%hdkh = 0._dp
928 :
929 504 : m = MAXVAL(basis%nprim)
930 792 : ALLOCATE (tp(m, m, 0:lmat), sp(m, m, 0:lmat), vp(m, m, 0:lmat), pvp(m, m, 0:lmat))
931 334608 : tp = 0._dp
932 334608 : sp = 0._dp
933 334608 : vp = 0._dp
934 334608 : pvp = 0._dp
935 :
936 72 : SELECT CASE (basis%basis_type)
937 : CASE DEFAULT
938 0 : CPABORT("")
939 : CASE (GTO_BASIS, CGTO_BASIS)
940 :
941 504 : DO l = 0, lmat
942 432 : m = basis%nprim(l)
943 504 : IF (m > 0) THEN
944 272 : CALL sg_kinetic(tp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
945 272 : CALL sg_overlap(sp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
946 272 : IF (PRESENT(alpha)) THEN
947 36 : CALL sg_erfc(vp(1:m, 1:m, l), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
948 : ELSE
949 236 : CALL sg_nuclear(vp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
950 : END IF
951 272 : CALL sg_kinnuc(pvp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
952 223656 : vp(1:m, 1:m, l) = -zcore*vp(1:m, 1:m, l)
953 223656 : pvp(1:m, 1:m, l) = -zcore*pvp(1:m, 1:m, l)
954 : END IF
955 : END DO
956 :
957 : CASE (STO_BASIS)
958 0 : CPABORT("")
959 : CASE (NUM_BASIS)
960 72 : CPABORT("")
961 : END SELECT
962 :
963 72 : CALL dkh_integrals(integrals, basis, dkhorder, sp, tp, vp, pvp)
964 :
965 72 : integrals%dkhstat = 1
966 :
967 72 : DEALLOCATE (tp, sp, vp, pvp)
968 :
969 : ELSE
970 0 : CPASSERT(ASSOCIATED(integrals%hdkh))
971 : END IF
972 :
973 : END SELECT
974 :
975 9908 : CALL timestop(handle)
976 :
977 9908 : END SUBROUTINE atom_relint_setup
978 :
979 : ! **************************************************************************************************
980 : !> \brief ...
981 : !> \param integrals ...
982 : !> \param basis ...
983 : !> \param order ...
984 : !> \param sp ...
985 : !> \param tp ...
986 : !> \param vp ...
987 : !> \param pvp ...
988 : ! **************************************************************************************************
989 72 : SUBROUTINE dkh_integrals(integrals, basis, order, sp, tp, vp, pvp)
990 : TYPE(atom_integrals), INTENT(INOUT) :: integrals
991 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
992 : INTEGER, INTENT(IN) :: order
993 : REAL(dp), DIMENSION(:, :, 0:) :: sp, tp, vp, pvp
994 :
995 : INTEGER :: l, m, n
996 72 : REAL(dp), DIMENSION(:, :, :), POINTER :: hdkh
997 :
998 0 : CPASSERT(order >= 0)
999 :
1000 72 : hdkh => integrals%hdkh
1001 :
1002 504 : DO l = 0, lmat
1003 432 : n = integrals%n(l)
1004 432 : m = basis%nprim(l)
1005 504 : IF (n > 0) THEN
1006 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)
1007 272 : SELECT CASE (basis%basis_type)
1008 : CASE DEFAULT
1009 0 : CPABORT("")
1010 : CASE (GTO_BASIS)
1011 222 : CPASSERT(n == m)
1012 219890 : integrals%hdkh(1:n, 1:n, l) = tp(1:n, 1:n, l) + vp(1:n, 1:n, l)
1013 : CASE (CGTO_BASIS)
1014 50 : CALL contract2(integrals%hdkh(1:n, 1:n, l), tp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1015 50 : CALL contract2add(integrals%hdkh(1:n, 1:n, l), vp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1016 : CASE (STO_BASIS)
1017 0 : CPABORT("")
1018 : CASE (NUM_BASIS)
1019 272 : CPABORT("")
1020 : END SELECT
1021 : ELSE
1022 160 : integrals%hdkh(1:n, 1:n, l) = 0._dp
1023 : END IF
1024 : END DO
1025 :
1026 72 : END SUBROUTINE dkh_integrals
1027 :
1028 : ! **************************************************************************************************
1029 : !> \brief Calculate overlap matrix between two atomic basis sets.
1030 : !> \param ovlap overlap matrix <chi_{a,l} | chi_{b,l}>
1031 : !> \param basis_a first basis set
1032 : !> \param basis_b second basis set
1033 : ! **************************************************************************************************
1034 1 : SUBROUTINE atom_basis_projection_overlap(ovlap, basis_a, basis_b)
1035 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: ovlap
1036 : TYPE(atom_basis_type), INTENT(IN) :: basis_a, basis_b
1037 :
1038 : INTEGER :: i, j, l, ma, mb, na, nb
1039 : LOGICAL :: ebas
1040 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
1041 :
1042 1 : ebas = (basis_a%basis_type == basis_b%basis_type)
1043 :
1044 127 : ovlap = 0.0_dp
1045 :
1046 1 : IF (ebas) THEN
1047 0 : SELECT CASE (basis_a%basis_type)
1048 : CASE DEFAULT
1049 0 : CPABORT("")
1050 : CASE (GTO_BASIS)
1051 0 : DO l = 0, lmat
1052 0 : na = basis_a%nbas(l)
1053 0 : nb = basis_b%nbas(l)
1054 0 : CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l))
1055 : END DO
1056 : CASE (CGTO_BASIS)
1057 7 : DO l = 0, lmat
1058 6 : na = basis_a%nbas(l)
1059 6 : nb = basis_b%nbas(l)
1060 6 : ma = basis_a%nprim(l)
1061 6 : mb = basis_b%nprim(l)
1062 18 : ALLOCATE (omat(ma, mb))
1063 6 : CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l))
1064 6 : ovlap(1:na, 1:nb, l) = MATMUL(TRANSPOSE(basis_a%cm(1:ma, 1:na, l)), &
1065 1277 : MATMUL(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l)))
1066 7 : DEALLOCATE (omat)
1067 : END DO
1068 : CASE (STO_BASIS)
1069 0 : DO l = 0, lmat
1070 0 : na = basis_a%nbas(l)
1071 0 : nb = basis_b%nbas(l)
1072 : CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), &
1073 0 : basis_a%ns(1:na, l), basis_b%as(1:nb, l))
1074 : END DO
1075 : CASE (NUM_BASIS)
1076 0 : CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
1077 1 : DO l = 0, lmat
1078 0 : na = basis_a%nbas(l)
1079 0 : nb = basis_b%nbas(l)
1080 0 : DO j = 1, nb
1081 0 : DO i = 1, na
1082 0 : ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1083 : END DO
1084 : END DO
1085 : END DO
1086 : END SELECT
1087 : ELSE
1088 0 : CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
1089 0 : DO l = 0, lmat
1090 0 : na = basis_a%nbas(l)
1091 0 : nb = basis_b%nbas(l)
1092 0 : DO j = 1, nb
1093 0 : DO i = 1, na
1094 0 : ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1095 : END DO
1096 : END DO
1097 : END DO
1098 : END IF
1099 :
1100 1 : END SUBROUTINE atom_basis_projection_overlap
1101 :
1102 : ! **************************************************************************************************
1103 : !> \brief Release memory allocated for atomic integrals (valence electrons).
1104 : !> \param integrals atomic integrals
1105 : ! **************************************************************************************************
1106 11150 : SUBROUTINE atom_int_release(integrals)
1107 : TYPE(atom_integrals), INTENT(INOUT) :: integrals
1108 :
1109 : INTEGER :: ll
1110 :
1111 11150 : IF (ASSOCIATED(integrals%ovlp)) THEN
1112 11150 : DEALLOCATE (integrals%ovlp)
1113 : END IF
1114 11150 : IF (ASSOCIATED(integrals%kin)) THEN
1115 11150 : DEALLOCATE (integrals%kin)
1116 : END IF
1117 11150 : IF (ASSOCIATED(integrals%conf)) THEN
1118 8584 : DEALLOCATE (integrals%conf)
1119 : END IF
1120 1126150 : DO ll = 1, SIZE(integrals%ceri)
1121 1115000 : IF (ASSOCIATED(integrals%ceri(ll)%int)) THEN
1122 1218 : DEALLOCATE (integrals%ceri(ll)%int)
1123 : END IF
1124 1126150 : IF (ASSOCIATED(integrals%eeri(ll)%int)) THEN
1125 2128 : DEALLOCATE (integrals%eeri(ll)%int)
1126 : END IF
1127 : END DO
1128 11150 : IF (ASSOCIATED(integrals%utrans)) THEN
1129 11150 : DEALLOCATE (integrals%utrans)
1130 : END IF
1131 11150 : IF (ASSOCIATED(integrals%uptrans)) THEN
1132 11150 : DEALLOCATE (integrals%uptrans)
1133 : END IF
1134 :
1135 11150 : integrals%status = 0
1136 :
1137 11150 : END SUBROUTINE atom_int_release
1138 :
1139 : ! **************************************************************************************************
1140 : !> \brief Release memory allocated for atomic integrals (core electrons).
1141 : !> \param integrals atomic integrals
1142 : ! **************************************************************************************************
1143 11480 : SUBROUTINE atom_ppint_release(integrals)
1144 : TYPE(atom_integrals), INTENT(INOUT) :: integrals
1145 :
1146 11480 : IF (ASSOCIATED(integrals%hnl)) THEN
1147 11480 : DEALLOCATE (integrals%hnl)
1148 : END IF
1149 11480 : IF (ASSOCIATED(integrals%core)) THEN
1150 11480 : DEALLOCATE (integrals%core)
1151 : END IF
1152 11480 : IF (ASSOCIATED(integrals%clsd)) THEN
1153 11480 : DEALLOCATE (integrals%clsd)
1154 : END IF
1155 :
1156 11480 : integrals%ppstat = 0
1157 :
1158 11480 : END SUBROUTINE atom_ppint_release
1159 :
1160 : ! **************************************************************************************************
1161 : !> \brief Release memory allocated for atomic integrals (relativistic effects).
1162 : !> \param integrals atomic integrals
1163 : ! **************************************************************************************************
1164 11148 : SUBROUTINE atom_relint_release(integrals)
1165 : TYPE(atom_integrals), INTENT(INOUT) :: integrals
1166 :
1167 11148 : IF (ASSOCIATED(integrals%tzora)) THEN
1168 28 : DEALLOCATE (integrals%tzora)
1169 : END IF
1170 11148 : IF (ASSOCIATED(integrals%hdkh)) THEN
1171 72 : DEALLOCATE (integrals%hdkh)
1172 : END IF
1173 :
1174 11148 : integrals%zorastat = 0
1175 11148 : integrals%dkhstat = 0
1176 :
1177 11148 : END SUBROUTINE atom_relint_release
1178 :
1179 : ! **************************************************************************************************
1180 : !> \brief Calculate model potential. It is useful to guess atomic orbitals.
1181 : !> \param modpot model potential
1182 : !> \param grid atomic radial grid
1183 : !> \param zcore nuclear charge
1184 : ! **************************************************************************************************
1185 72 : SUBROUTINE calculate_model_potential(modpot, grid, zcore)
1186 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: modpot
1187 : TYPE(grid_atom_type), INTENT(IN) :: grid
1188 : REAL(dp), INTENT(IN) :: zcore
1189 :
1190 : INTEGER :: i, ii, l, ll, n, z
1191 72 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pot, rho
1192 : TYPE(atom_state) :: state
1193 :
1194 72 : n = SIZE(modpot)
1195 288 : ALLOCATE (rho(n), pot(n))
1196 :
1197 : ! fill default occupation
1198 5112 : state%core = 0._dp
1199 5112 : state%occ = 0._dp
1200 72 : state%multiplicity = -1
1201 72 : z = NINT(zcore)
1202 360 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
1203 360 : IF (ptable(z)%e_conv(l) /= 0) THEN
1204 156 : state%maxl_occ = l
1205 156 : ll = 2*(2*l + 1)
1206 360 : DO i = 1, SIZE(state%occ, 2)
1207 360 : ii = ptable(z)%e_conv(l) - (i - 1)*ll
1208 360 : IF (ii <= ll) THEN
1209 156 : state%occ(l, i) = ii
1210 156 : EXIT
1211 : ELSE
1212 204 : state%occ(l, i) = ll
1213 : END IF
1214 : END DO
1215 : END IF
1216 : END DO
1217 :
1218 13410 : modpot = -zcore/grid%rad(:)
1219 :
1220 : ! Coulomb potential
1221 72 : CALL slater_density(rho, pot, NINT(zcore), state, grid)
1222 72 : CALL coulomb_potential_numeric(pot, rho, grid)
1223 13410 : modpot = modpot + pot
1224 :
1225 : ! XC potential
1226 72 : CALL wigner_slater_functional(rho, pot)
1227 13410 : modpot = modpot + pot
1228 :
1229 72 : DEALLOCATE (rho, pot)
1230 :
1231 26136 : END SUBROUTINE calculate_model_potential
1232 :
1233 14137 : END MODULE atom_operators
|
