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