1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief   Calculate the atomic operator matrices
8!> \author  jgh
9!> \date    03.03.2008
10!> \version 1.0
11!>
12! **************************************************************************************************
13MODULE atom_operators
14   USE ai_onecenter,                    ONLY: &
15        sg_coulomb, sg_erf, sg_erfc, sg_exchange, sg_gpot, sg_kinetic, sg_kinnuc, sg_nuclear, &
16        sg_overlap, sg_proj_ol, sto_kinetic, sto_nuclear, sto_overlap
17   USE atom_types,                      ONLY: &
18        atom_basis_gridrep, atom_basis_type, atom_compare_grids, atom_integrals, &
19        atom_potential_type, atom_state, cgto_basis, ecp_pseudo, gth_pseudo, gto_basis, lmat, &
20        no_pseudo, num_basis, release_atom_basis, sgp_pseudo, sto_basis, upf_pseudo
21   USE atom_utils,                      ONLY: &
22        atom_solve, contract2, contract2add, contract4, coulomb_potential_numeric, integrate_grid, &
23        numpot_matrix, slater_density, wigner_slater_functional
24   USE dkh_main,                        ONLY: dkh_atom_transformation
25   USE input_constants,                 ONLY: &
26        barrier_conf, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
27        do_sczoramp_atom, do_zoramp_atom, poly_conf
28   USE kinds,                           ONLY: dp
29   USE lapack,                          ONLY: lapack_sgesv
30   USE mathconstants,                   ONLY: gamma1,&
31                                              sqrt2
32   USE mathlib,                         ONLY: jacobi
33   USE periodic_table,                  ONLY: ptable
34   USE physcon,                         ONLY: c_light_au
35   USE qs_grid_atom,                    ONLY: grid_atom_type
36#include "./base/base_uses.f90"
37
38   IMPLICIT NONE
39
40   PRIVATE
41
42   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_operators'
43
44   PUBLIC :: atom_int_setup, atom_ppint_setup, atom_int_release, atom_ppint_release
45   PUBLIC :: atom_relint_setup, atom_relint_release, atom_basis_projection_overlap
46   PUBLIC :: calculate_model_potential
47
48CONTAINS
49
50! **************************************************************************************************
51!> \brief Set up atomic integrals.
52!> \param integrals     atomic integrals
53!> \param basis         atomic basis set
54!> \param potential     pseudo-potential
55!> \param eri_coulomb   setup one-centre Coulomb Electron Repulsion Integrals
56!> \param eri_exchange  setup one-centre exchange Electron Repulsion Integrals
57!> \param all_nu        compute integrals for all even integer parameters [0 .. 2*l]
58!>                      REDUNDANT, AS THIS SUBROUTINE IS NEVER INVOKED WITH all_nu = .TRUE.
59! **************************************************************************************************
60   SUBROUTINE atom_int_setup(integrals, basis, potential, &
61                             eri_coulomb, eri_exchange, all_nu)
62      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
63      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
64      TYPE(atom_potential_type), INTENT(IN), OPTIONAL    :: potential
65      LOGICAL, INTENT(IN), OPTIONAL                      :: eri_coulomb, eri_exchange, all_nu
66
67      CHARACTER(len=*), PARAMETER :: routineN = 'atom_int_setup', routineP = moduleN//':'//routineN
68
69      INTEGER                                            :: handle, i, ii, info, ipiv(1000), l, l1, &
70                                                            l2, ll, lwork, m, m1, m2, mm1, mm2, n, &
71                                                            n1, n2, nn1, nn2, nu, nx
72      REAL(KIND=dp)                                      :: om, rc, ron, sc, x
73      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cpot, w, work
74      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: omat, vmat
75      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: eri
76
77      CALL timeset(routineN, handle)
78
79      IF (integrals%status == 0) THEN
80         n = MAXVAL(basis%nbas)
81         integrals%n = basis%nbas
82
83         IF (PRESENT(eri_coulomb)) THEN
84            integrals%eri_coulomb = eri_coulomb
85         ELSE
86            integrals%eri_coulomb = .FALSE.
87         END IF
88         IF (PRESENT(eri_exchange)) THEN
89            integrals%eri_exchange = eri_exchange
90         ELSE
91            integrals%eri_exchange = .FALSE.
92         END IF
93         IF (PRESENT(all_nu)) THEN
94            integrals%all_nu = all_nu
95         ELSE
96            integrals%all_nu = .FALSE.
97         END IF
98
99         NULLIFY (integrals%ovlp, integrals%kin, integrals%core, integrals%conf)
100         DO ll = 1, SIZE(integrals%ceri)
101            NULLIFY (integrals%ceri(ll)%int, integrals%eeri(ll)%int)
102         END DO
103
104         ALLOCATE (integrals%ovlp(n, n, 0:lmat))
105         integrals%ovlp = 0._dp
106
107         ALLOCATE (integrals%kin(n, n, 0:lmat))
108         integrals%kin = 0._dp
109
110         integrals%status = 1
111
112         IF (PRESENT(potential)) THEN
113            IF (potential%confinement) THEN
114               ALLOCATE (integrals%conf(n, n, 0:lmat))
115               integrals%conf = 0._dp
116               m = basis%grid%nr
117               ALLOCATE (cpot(1:m))
118               IF (potential%conf_type == poly_conf) THEN
119                  rc = potential%rcon
120                  sc = potential%scon
121                  cpot(1:m) = (basis%grid%rad(1:m)/rc)**sc
122               ELSEIF (potential%conf_type == barrier_conf) THEN
123                  om = potential%rcon
124                  ron = potential%scon
125                  rc = ron + om
126                  DO i = 1, m
127                     IF (basis%grid%rad(i) < ron) THEN
128                        cpot(i) = 0.0_dp
129                     ELSEIF (basis%grid%rad(i) < rc) THEN
130                        x = (basis%grid%rad(i) - ron)/om
131                        x = 1._dp - x
132                        cpot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
133                        x = (rc - basis%grid%rad(i))**2/om/(basis%grid%rad(i) - ron)
134                        cpot(i) = cpot(i)*x
135                     ELSE
136                        cpot(i) = 1.0_dp
137                     END IF
138                  END DO
139               ELSE
140                  CPABORT("")
141               END IF
142               CALL numpot_matrix(integrals%conf, cpot, basis, 0)
143               DEALLOCATE (cpot)
144            END IF
145         END IF
146
147         SELECT CASE (basis%basis_type)
148         CASE DEFAULT
149            CPABORT("")
150         CASE (GTO_BASIS)
151            DO l = 0, lmat
152               n = integrals%n(l)
153               CALL sg_overlap(integrals%ovlp(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
154               CALL sg_kinetic(integrals%kin(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
155            END DO
156            IF (integrals%eri_coulomb) THEN
157               ll = 0
158               DO l1 = 0, lmat
159                  n1 = integrals%n(l1)
160                  nn1 = (n1*(n1 + 1))/2
161                  DO l2 = 0, l1
162                     n2 = integrals%n(l2)
163                     nn2 = (n2*(n2 + 1))/2
164                     IF (integrals%all_nu) THEN
165                        nx = MIN(2*l1, 2*l2)
166                     ELSE
167                        nx = 0
168                     END IF
169                     DO nu = 0, nx, 2
170                        ll = ll + 1
171                        CPASSERT(ll <= SIZE(integrals%ceri))
172                        ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
173                        integrals%ceri(ll)%int = 0._dp
174                        eri => integrals%ceri(ll)%int
175                        CALL sg_coulomb(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
176                     END DO
177                  END DO
178               END DO
179            END IF
180            IF (integrals%eri_exchange) THEN
181               ll = 0
182               DO l1 = 0, lmat
183                  n1 = integrals%n(l1)
184                  nn1 = (n1*(n1 + 1))/2
185                  DO l2 = 0, l1
186                     n2 = integrals%n(l2)
187                     nn2 = (n2*(n2 + 1))/2
188                     DO nu = ABS(l1 - l2), l1 + l2, 2
189                        ll = ll + 1
190                        CPASSERT(ll <= SIZE(integrals%eeri))
191                        ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
192                        integrals%eeri(ll)%int = 0._dp
193                        eri => integrals%eeri(ll)%int
194                        CALL sg_exchange(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
195                     END DO
196                  END DO
197               END DO
198            END IF
199         CASE (CGTO_BASIS)
200            DO l = 0, lmat
201               n = integrals%n(l)
202               m = basis%nprim(l)
203               IF (n > 0 .AND. m > 0) THEN
204                  ALLOCATE (omat(m, m))
205
206                  CALL sg_overlap(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
207                  CALL contract2(integrals%ovlp(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
208                  CALL sg_kinetic(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
209                  CALL contract2(integrals%kin(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
210                  DEALLOCATE (omat)
211               END IF
212            END DO
213            IF (integrals%eri_coulomb) THEN
214               ll = 0
215               DO l1 = 0, lmat
216                  n1 = integrals%n(l1)
217                  nn1 = (n1*(n1 + 1))/2
218                  m1 = basis%nprim(l1)
219                  mm1 = (m1*(m1 + 1))/2
220                  DO l2 = 0, l1
221                     n2 = integrals%n(l2)
222                     nn2 = (n2*(n2 + 1))/2
223                     m2 = basis%nprim(l2)
224                     mm2 = (m2*(m2 + 1))/2
225                     IF (integrals%all_nu) THEN
226                        nx = MIN(2*l1, 2*l2)
227                     ELSE
228                        nx = 0
229                     END IF
230                     DO nu = 0, nx, 2
231                        ll = ll + 1
232                        CPASSERT(ll <= SIZE(integrals%ceri))
233                        ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
234                        integrals%ceri(ll)%int = 0._dp
235                        ALLOCATE (omat(mm1, mm2))
236
237                        eri => integrals%ceri(ll)%int
238                        CALL sg_coulomb(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
239                        CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))
240
241                        DEALLOCATE (omat)
242                     END DO
243                  END DO
244               END DO
245            END IF
246            IF (integrals%eri_exchange) THEN
247               ll = 0
248               DO l1 = 0, lmat
249                  n1 = integrals%n(l1)
250                  nn1 = (n1*(n1 + 1))/2
251                  m1 = basis%nprim(l1)
252                  mm1 = (m1*(m1 + 1))/2
253                  DO l2 = 0, l1
254                     n2 = integrals%n(l2)
255                     nn2 = (n2*(n2 + 1))/2
256                     m2 = basis%nprim(l2)
257                     mm2 = (m2*(m2 + 1))/2
258                     DO nu = ABS(l1 - l2), l1 + l2, 2
259                        ll = ll + 1
260                        CPASSERT(ll <= SIZE(integrals%eeri))
261                        ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
262                        integrals%eeri(ll)%int = 0._dp
263                        ALLOCATE (omat(mm1, mm2))
264
265                        eri => integrals%eeri(ll)%int
266                        CALL sg_exchange(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
267                        CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))
268
269                        DEALLOCATE (omat)
270                     END DO
271                  END DO
272               END DO
273            END IF
274         CASE (STO_BASIS)
275            DO l = 0, lmat
276               n = integrals%n(l)
277               CALL sto_overlap(integrals%ovlp(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
278                                basis%ns(1:n, l), basis%as(1:n, l))
279               CALL sto_kinetic(integrals%kin(1:n, 1:n, l), l, basis%ns(1:n, l), basis%as(1:n, l), &
280                                basis%ns(1:n, l), basis%as(1:n, l))
281            END DO
282            CPASSERT(.NOT. integrals%eri_coulomb)
283            CPASSERT(.NOT. integrals%eri_exchange)
284         CASE (NUM_BASIS)
285            CPABORT("")
286         END SELECT
287
288         ! setup transformation matrix to get an orthogonal basis, remove linear dependencies
289         NULLIFY (integrals%utrans, integrals%uptrans)
290         n = MAXVAL(basis%nbas)
291         ALLOCATE (integrals%utrans(n, n, 0:lmat), integrals%uptrans(n, n, 0:lmat))
292         integrals%utrans = 0._dp
293         integrals%uptrans = 0._dp
294         integrals%nne = integrals%n
295         lwork = 10*n
296         ALLOCATE (omat(n, n), w(n), vmat(n, n), work(lwork))
297         DO l = 0, lmat
298            n = integrals%n(l)
299            IF (n > 0) THEN
300               omat(1:n, 1:n) = integrals%ovlp(1:n, 1:n, l)
301               CALL jacobi(omat(1:n, 1:n), w(1:n), vmat(1:n, 1:n))
302               omat(1:n, 1:n) = vmat(1:n, 1:n)
303               ii = 0
304               DO i = 1, n
305                  IF (w(i) > basis%eps_eig) THEN
306                     ii = ii + 1
307                     integrals%utrans(1:n, ii, l) = omat(1:n, i)/SQRT(w(i))
308                  END IF
309               END DO
310               integrals%nne(l) = ii
311               IF (ii > 0) THEN
312                  omat(1:ii, 1:ii) = MATMUL(TRANSPOSE(integrals%utrans(1:n, 1:ii, l)), integrals%utrans(1:n, 1:ii, l))
313                  DO i = 1, ii
314                     integrals%uptrans(i, i, l) = 1._dp
315                  ENDDO
316                  CALL lapack_sgesv(ii, ii, omat(1:ii, 1:ii), ii, ipiv, integrals%uptrans(1:ii, 1:ii, l), ii, info)
317                  CPASSERT(info == 0)
318               END IF
319            END IF
320         END DO
321         DEALLOCATE (omat, vmat, w, work)
322
323      END IF
324
325      CALL timestop(handle)
326
327   END SUBROUTINE atom_int_setup
328
329! **************************************************************************************************
330!> \brief ...
331!> \param integrals ...
332!> \param basis ...
333!> \param potential ...
334! **************************************************************************************************
335   SUBROUTINE atom_ppint_setup(integrals, basis, potential)
336      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
337      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
338      TYPE(atom_potential_type), INTENT(IN)              :: potential
339
340      CHARACTER(len=*), PARAMETER :: routineN = 'atom_ppint_setup', &
341         routineP = moduleN//':'//routineN
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=*), PARAMETER :: routineN = 'upfint_setup', routineP = moduleN//':'//routineN
652
653      CHARACTER(len=4)                                   :: ptype
654      INTEGER                                            :: i, j, k1, k2, la, lb, m, n
655      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: spot
656      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: spmat
657      TYPE(atom_basis_type)                              :: gbasis
658
659      ! get basis representation on UPF grid
660      CALL atom_basis_gridrep(basis, gbasis, potential%upf_pot%r, potential%upf_pot%rab)
661
662      ! local pseudopotential
663      integrals%core = 0._dp
664      CALL numpot_matrix(integrals%core, potential%upf_pot%vlocal, gbasis, 0)
665
666      ptype = ADJUSTL(TRIM(potential%upf_pot%pseudo_type))
667      IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
668         ! non local pseudopotential
669         n = MAXVAL(integrals%n(:))
670         m = potential%upf_pot%number_of_proj
671         ALLOCATE (spmat(n, m))
672         spmat = 0.0_dp
673         DO i = 1, m
674            la = potential%upf_pot%lbeta(i)
675            DO j = 1, gbasis%nbas(la)
676               spmat(j, i) = integrate_grid(potential%upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
677            END DO
678         END DO
679         DO i = 1, m
680            la = potential%upf_pot%lbeta(i)
681            DO j = 1, m
682               lb = potential%upf_pot%lbeta(j)
683               IF (la == lb) THEN
684                  DO k1 = 1, gbasis%nbas(la)
685                     DO k2 = 1, gbasis%nbas(la)
686                        integrals%hnl(k1, k2, la) = integrals%hnl(k1, k2, la) + &
687                                                    spmat(k1, i)*potential%upf_pot%dion(i, j)*spmat(k2, j)
688                     END DO
689                  END DO
690               END IF
691            END DO
692         END DO
693         DEALLOCATE (spmat)
694      ELSE IF (ptype(1:2) == "SL") THEN
695         ! semi local pseudopotential
696         DO la = 0, potential%upf_pot%l_max
697            IF (la == potential%upf_pot%l_local) CYCLE
698            m = SIZE(potential%upf_pot%vsemi(:, la + 1))
699            ALLOCATE (spot(m))
700            spot(:) = potential%upf_pot%vsemi(:, la + 1) - potential%upf_pot%vlocal(:)
701            n = basis%nbas(la)
702            DO i = 1, n
703               DO j = i, n
704                  integrals%core(i, j, la) = integrals%core(i, j, la) + &
705                                             integrate_grid(spot(:), &
706                                                            gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
707                  integrals%core(j, i, la) = integrals%core(i, j, la)
708               END DO
709            END DO
710            DEALLOCATE (spot)
711         END DO
712      ELSE
713         CPABORT("Pseudopotential type: ["//ADJUSTL(TRIM(ptype))//"] not known")
714      END IF
715
716      ! release basis representation on UPF grid
717      CALL release_atom_basis(gbasis)
718
719   END SUBROUTINE upfint_setup
720
721! **************************************************************************************************
722!> \brief ...
723!> \param integrals ...
724!> \param basis ...
725!> \param potential ...
726! **************************************************************************************************
727   SUBROUTINE sgpint_setup(integrals, basis, potential)
728      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
729      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
730      TYPE(atom_potential_type), INTENT(IN)              :: potential
731
732      CHARACTER(len=*), PARAMETER :: routineN = 'sgpint_setup', routineP = moduleN//':'//routineN
733
734      INTEGER                                            :: i, ia, j, l, m, n, na
735      REAL(KIND=dp)                                      :: a, c, rc, zval
736      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cpot, pgauss
737      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: qmat
738      REAL(KIND=dp), DIMENSION(:), POINTER               :: rad
739
740      rad => basis%grid%rad
741      m = basis%grid%nr
742
743      ! local pseudopotential
744      integrals%core = 0._dp
745      ALLOCATE (cpot(m))
746      cpot = 0.0_dp
747      zval = potential%sgp_pot%zion
748      DO i = 1, m
749         rc = rad(i)/potential%sgp_pot%ac_local/SQRT(2.0_dp)
750         cpot(i) = cpot(i) - zval/rad(i)*erf(rc)
751      END DO
752      DO i = 1, potential%sgp_pot%n_local
753         cpot(:) = cpot(:) + potential%sgp_pot%c_local(i)*EXP(-potential%sgp_pot%a_local(i)*rad(:)**2)
754      END DO
755      CALL numpot_matrix(integrals%core, cpot, basis, 0)
756      DEALLOCATE (cpot)
757
758      ! nonlocal pseudopotential
759      integrals%hnl = 0.0_dp
760      IF (potential%sgp_pot%has_nonlocal) THEN
761         ALLOCATE (pgauss(1:m))
762         n = potential%sgp_pot%n_nonlocal
763         !
764         DO l = 0, potential%sgp_pot%lmax
765            CPASSERT(l <= UBOUND(basis%nbas, 1))
766            IF (.NOT. potential%sgp_pot%is_nonlocal(l)) CYCLE
767            ! overlap (a|p)
768            na = basis%nbas(l)
769            ALLOCATE (qmat(na, n))
770            DO i = 1, n
771               pgauss(:) = 0.0_dp
772               DO j = 1, n
773                  a = potential%sgp_pot%a_nonlocal(j)
774                  c = potential%sgp_pot%c_nonlocal(j, i, l)
775                  pgauss(:) = pgauss(:) + c*EXP(-a*rad(:)**2)*rad(:)**l
776               END DO
777               DO ia = 1, na
778                  qmat(ia, i) = SUM(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:))
779               END DO
780            END DO
781            DO i = 1, na
782               DO j = i, na
783                  DO ia = 1, n
784                     integrals%hnl(i, j, l) = integrals%hnl(i, j, l) &
785                                              + qmat(i, ia)*qmat(j, ia)*potential%sgp_pot%h_nonlocal(ia, l)
786                  END DO
787                  integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
788               END DO
789            END DO
790            DEALLOCATE (qmat)
791         END DO
792         DEALLOCATE (pgauss)
793      END IF
794
795   END SUBROUTINE sgpint_setup
796
797! **************************************************************************************************
798!> \brief ...
799!> \param integrals ...
800!> \param basis ...
801!> \param reltyp ...
802!> \param zcore ...
803!> \param alpha ...
804! **************************************************************************************************
805   SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
806      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
807      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
808      INTEGER, INTENT(IN)                                :: reltyp
809      REAL(dp), INTENT(IN)                               :: zcore
810      REAL(dp), INTENT(IN), OPTIONAL                     :: alpha
811
812      CHARACTER(len=*), PARAMETER :: routineN = 'atom_relint_setup', &
813         routineP = moduleN//':'//routineN
814
815      INTEGER                                            :: dkhorder, handle, i, k1, k2, l, m, n, nl
816      REAL(dp)                                           :: ascal
817      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: cpot
818      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: modpot
819      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ener, sps
820      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: hmat, pvp, sp, tp, vp, wfn
821
822      CALL timeset(routineN, handle)
823
824      SELECT CASE (reltyp)
825      CASE DEFAULT
826         CPABORT("")
827      CASE (do_nonrel_atom, do_zoramp_atom, do_sczoramp_atom)
828         dkhorder = -1
829      CASE (do_dkh0_atom)
830         dkhorder = 0
831      CASE (do_dkh1_atom)
832         dkhorder = 1
833      CASE (do_dkh2_atom)
834         dkhorder = 2
835      CASE (do_dkh3_atom)
836         dkhorder = 3
837      END SELECT
838
839      SELECT CASE (reltyp)
840      CASE DEFAULT
841         CPABORT("")
842      CASE (do_nonrel_atom)
843         ! nothing to do
844         NULLIFY (integrals%tzora, integrals%hdkh)
845      CASE (do_zoramp_atom, do_sczoramp_atom)
846
847         NULLIFY (integrals%hdkh)
848
849         IF (integrals%zorastat == 0) THEN
850            n = MAXVAL(basis%nbas)
851            ALLOCATE (integrals%tzora(n, n, 0:lmat))
852            integrals%tzora = 0._dp
853            m = basis%grid%nr
854            ALLOCATE (modpot(1:m), cpot(1:m))
855            CALL calculate_model_potential(modpot, basis%grid, zcore)
856            ! Zora potential
857            cpot(1:m) = modpot(1:m)/(4._dp*c_light_au*c_light_au - 2._dp*modpot(1:m))
858            cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
859            CALL numpot_matrix(integrals%tzora, cpot, basis, 0)
860            DO l = 0, lmat
861               nl = basis%nbas(l)
862               integrals%tzora(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*integrals%tzora(1:nl, 1:nl, l)
863            END DO
864            cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
865            CALL numpot_matrix(integrals%tzora, cpot, basis, 2)
866            !
867            ! scaled ZORA
868            IF (reltyp == do_sczoramp_atom) THEN
869               ALLOCATE (hmat(n, n, 0:lmat), wfn(n, n, 0:lmat), ener(n, 0:lmat), pvp(n, n, 0:lmat), sps(n, n))
870               hmat(:, :, :) = integrals%kin + integrals%tzora
871               ! model potential
872               CALL numpot_matrix(hmat, modpot, basis, 0)
873               ! eigenvalues and eigenvectors
874               CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne, lmat)
875               ! relativistic kinetic energy
876               cpot(1:m) = c_light_au*c_light_au/(2._dp*c_light_au*c_light_au - modpot(1:m))**2
877               cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
878               pvp = 0.0_dp
879               CALL numpot_matrix(pvp, cpot, basis, 0)
880               DO l = 0, lmat
881                  nl = basis%nbas(l)
882                  pvp(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*pvp(1:nl, 1:nl, l)
883               END DO
884               cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
885               CALL numpot_matrix(pvp, cpot, basis, 2)
886               ! calculate psi*pvp*psi and the scaled orbital energies
887               ! actually, we directly calculate the energy difference
888               DO l = 0, lmat
889                  nl = basis%nbas(l)
890                  DO i = 1, integrals%nne(l)
891                     IF (ener(i, l) < 0._dp) THEN
892                        ascal = SUM(wfn(1:nl, i, l)*MATMUL(pvp(1:nl, 1:nl, l), wfn(1:nl, i, l)))
893                        ener(i, l) = ener(i, l)*ascal/(1.0_dp + ascal)
894                     ELSE
895                        ener(i, l) = 0.0_dp
896                     END IF
897                  END DO
898               END DO
899               ! correction term is calculated as a projector
900               hmat = 0.0_dp
901               DO l = 0, lmat
902                  nl = basis%nbas(l)
903                  DO i = 1, integrals%nne(l)
904                     DO k1 = 1, nl
905                        DO k2 = 1, nl
906                           hmat(k1, k2, l) = hmat(k1, k2, l) + ener(i, l)*wfn(k1, i, l)*wfn(k2, i, l)
907                        END DO
908                     END DO
909                  END DO
910                  ! transform with overlap matrix
911                  sps(1:nl, 1:nl) = MATMUL(integrals%ovlp(1:nl, 1:nl, l), &
912                                           MATMUL(hmat(1:nl, 1:nl, l), integrals%ovlp(1:nl, 1:nl, l)))
913                  ! add scaling correction to tzora
914                  integrals%tzora(1:nl, 1:nl, l) = integrals%tzora(1:nl, 1:nl, l) - sps(1:nl, 1:nl)
915               END DO
916
917               DEALLOCATE (hmat, wfn, ener, pvp, sps)
918            END IF
919            !
920            DEALLOCATE (modpot, cpot)
921
922            integrals%zorastat = 1
923
924         END IF
925
926      CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
927
928         NULLIFY (integrals%tzora)
929
930         IF (integrals%dkhstat == 0) THEN
931            n = MAXVAL(basis%nbas)
932            ALLOCATE (integrals%hdkh(n, n, 0:lmat))
933            integrals%hdkh = 0._dp
934
935            m = MAXVAL(basis%nprim)
936            ALLOCATE (tp(m, m, 0:lmat), sp(m, m, 0:lmat), vp(m, m, 0:lmat), pvp(m, m, 0:lmat))
937            tp = 0._dp
938            sp = 0._dp
939            vp = 0._dp
940            pvp = 0._dp
941
942            SELECT CASE (basis%basis_type)
943            CASE DEFAULT
944               CPABORT("")
945            CASE (GTO_BASIS, CGTO_BASIS)
946
947               DO l = 0, lmat
948                  m = basis%nprim(l)
949                  IF (m > 0) THEN
950                     CALL sg_kinetic(tp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
951                     CALL sg_overlap(sp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
952                     IF (PRESENT(alpha)) THEN
953                        CALL sg_erfc(vp(1:m, 1:m, l), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
954                     ELSE
955                        CALL sg_nuclear(vp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
956                     END IF
957                     CALL sg_kinnuc(pvp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
958                     vp(1:m, 1:m, l) = -zcore*vp(1:m, 1:m, l)
959                     pvp(1:m, 1:m, l) = -zcore*pvp(1:m, 1:m, l)
960                  END IF
961               END DO
962
963            CASE (STO_BASIS)
964               CPABORT("")
965            CASE (NUM_BASIS)
966               CPABORT("")
967            END SELECT
968
969            CALL dkh_integrals(integrals, basis, dkhorder, sp, tp, vp, pvp)
970
971            integrals%dkhstat = 1
972
973            DEALLOCATE (tp, sp, vp, pvp)
974
975         ELSE
976            CPASSERT(ASSOCIATED(integrals%hdkh))
977         END IF
978
979      END SELECT
980
981      CALL timestop(handle)
982
983   END SUBROUTINE atom_relint_setup
984
985! **************************************************************************************************
986!> \brief ...
987!> \param integrals ...
988!> \param basis ...
989!> \param order ...
990!> \param sp ...
991!> \param tp ...
992!> \param vp ...
993!> \param pvp ...
994! **************************************************************************************************
995   SUBROUTINE dkh_integrals(integrals, basis, order, sp, tp, vp, pvp)
996      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
997      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
998      INTEGER, INTENT(IN)                                :: order
999      REAL(dp), DIMENSION(:, :, 0:)                      :: sp, tp, vp, pvp
1000
1001      CHARACTER(len=*), PARAMETER :: routineN = 'dkh_integrals', routineP = moduleN//':'//routineN
1002
1003      INTEGER                                            :: l, m, n
1004      REAL(dp), DIMENSION(:, :, :), POINTER              :: hdkh
1005
1006      CPASSERT(order >= 0)
1007
1008      hdkh => integrals%hdkh
1009
1010      DO l = 0, lmat
1011         n = integrals%n(l)
1012         m = basis%nprim(l)
1013         IF (n > 0) THEN
1014            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)
1015            SELECT CASE (basis%basis_type)
1016            CASE DEFAULT
1017               CPABORT("")
1018            CASE (GTO_BASIS)
1019               CPASSERT(n == m)
1020               integrals%hdkh(1:n, 1:n, l) = tp(1:n, 1:n, l) + vp(1:n, 1:n, l)
1021            CASE (CGTO_BASIS)
1022               CALL contract2(integrals%hdkh(1:n, 1:n, l), tp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1023               CALL contract2add(integrals%hdkh(1:n, 1:n, l), vp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
1024            CASE (STO_BASIS)
1025               CPABORT("")
1026            CASE (NUM_BASIS)
1027               CPABORT("")
1028            END SELECT
1029         ELSE
1030            integrals%hdkh(1:n, 1:n, l) = 0._dp
1031         END IF
1032      END DO
1033
1034   END SUBROUTINE dkh_integrals
1035
1036! **************************************************************************************************
1037!> \brief Calculate overlap matrix between two atomic basis sets.
1038!> \param ovlap    overlap matrix <chi_{a,l} | chi_{b,l}>
1039!> \param basis_a  first basis set
1040!> \param basis_b  second basis set
1041! **************************************************************************************************
1042   SUBROUTINE atom_basis_projection_overlap(ovlap, basis_a, basis_b)
1043      REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(OUT)    :: ovlap
1044      TYPE(atom_basis_type), INTENT(IN)                  :: basis_a, basis_b
1045
1046      INTEGER                                            :: i, j, l, ma, mb, na, nb
1047      LOGICAL                                            :: ebas
1048      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: omat
1049
1050      ebas = (basis_a%basis_type == basis_b%basis_type)
1051
1052      ovlap = 0.0_dp
1053
1054      IF (ebas) THEN
1055         SELECT CASE (basis_a%basis_type)
1056         CASE DEFAULT
1057            CPABORT("")
1058         CASE (GTO_BASIS)
1059            DO l = 0, lmat
1060               na = basis_a%nbas(l)
1061               nb = basis_b%nbas(l)
1062               CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l))
1063            END DO
1064         CASE (CGTO_BASIS)
1065            DO l = 0, lmat
1066               na = basis_a%nbas(l)
1067               nb = basis_b%nbas(l)
1068               ma = basis_a%nprim(l)
1069               mb = basis_b%nprim(l)
1070               ALLOCATE (omat(ma, mb))
1071               CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l))
1072               ovlap(1:na, 1:nb, l) = MATMUL(TRANSPOSE(basis_a%cm(1:ma, 1:na, l)), &
1073                                             MATMUL(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l)))
1074               DEALLOCATE (omat)
1075            END DO
1076         CASE (STO_BASIS)
1077            DO l = 0, lmat
1078               na = basis_a%nbas(l)
1079               nb = basis_b%nbas(l)
1080               CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), &
1081                                basis_a%ns(1:na, l), basis_b%as(1:nb, l))
1082            END DO
1083         CASE (NUM_BASIS)
1084            CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
1085            DO l = 0, lmat
1086               na = basis_a%nbas(l)
1087               nb = basis_b%nbas(l)
1088               DO j = 1, nb
1089                  DO i = 1, na
1090                     ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1091                  END DO
1092               END DO
1093            END DO
1094         END SELECT
1095      ELSE
1096         CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
1097         DO l = 0, lmat
1098            na = basis_a%nbas(l)
1099            nb = basis_b%nbas(l)
1100            DO j = 1, nb
1101               DO i = 1, na
1102                  ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
1103               END DO
1104            END DO
1105         END DO
1106      END IF
1107
1108   END SUBROUTINE atom_basis_projection_overlap
1109
1110! **************************************************************************************************
1111!> \brief Release memory allocated for atomic integrals (valence electrons).
1112!> \param integrals  atomic integrals
1113! **************************************************************************************************
1114   SUBROUTINE atom_int_release(integrals)
1115      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
1116
1117      CHARACTER(len=*), PARAMETER :: routineN = 'atom_int_release', &
1118         routineP = moduleN//':'//routineN
1119
1120      INTEGER                                            :: ll
1121
1122      IF (ASSOCIATED(integrals%ovlp)) THEN
1123         DEALLOCATE (integrals%ovlp)
1124      END IF
1125      IF (ASSOCIATED(integrals%kin)) THEN
1126         DEALLOCATE (integrals%kin)
1127      END IF
1128      IF (ASSOCIATED(integrals%conf)) THEN
1129         DEALLOCATE (integrals%conf)
1130      END IF
1131      DO ll = 1, SIZE(integrals%ceri)
1132         IF (ASSOCIATED(integrals%ceri(ll)%int)) THEN
1133            DEALLOCATE (integrals%ceri(ll)%int)
1134         END IF
1135         IF (ASSOCIATED(integrals%eeri(ll)%int)) THEN
1136            DEALLOCATE (integrals%eeri(ll)%int)
1137         END IF
1138      END DO
1139      IF (ASSOCIATED(integrals%utrans)) THEN
1140         DEALLOCATE (integrals%utrans)
1141      END IF
1142      IF (ASSOCIATED(integrals%uptrans)) THEN
1143         DEALLOCATE (integrals%uptrans)
1144      END IF
1145
1146      integrals%status = 0
1147
1148   END SUBROUTINE atom_int_release
1149
1150! **************************************************************************************************
1151!> \brief Release memory allocated for atomic integrals (core electrons).
1152!> \param integrals  atomic integrals
1153! **************************************************************************************************
1154   SUBROUTINE atom_ppint_release(integrals)
1155      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
1156
1157      CHARACTER(len=*), PARAMETER :: routineN = 'atom_ppint_release', &
1158         routineP = moduleN//':'//routineN
1159
1160      IF (ASSOCIATED(integrals%hnl)) THEN
1161         DEALLOCATE (integrals%hnl)
1162      END IF
1163      IF (ASSOCIATED(integrals%core)) THEN
1164         DEALLOCATE (integrals%core)
1165      END IF
1166      IF (ASSOCIATED(integrals%clsd)) THEN
1167         DEALLOCATE (integrals%clsd)
1168      END IF
1169
1170      integrals%ppstat = 0
1171
1172   END SUBROUTINE atom_ppint_release
1173
1174! **************************************************************************************************
1175!> \brief  Release memory allocated for atomic integrals (relativistic effects).
1176!> \param integrals atomic integrals
1177! **************************************************************************************************
1178   SUBROUTINE atom_relint_release(integrals)
1179      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
1180
1181      CHARACTER(len=*), PARAMETER :: routineN = 'atom_relint_release', &
1182         routineP = moduleN//':'//routineN
1183
1184      IF (ASSOCIATED(integrals%tzora)) THEN
1185         DEALLOCATE (integrals%tzora)
1186      END IF
1187      IF (ASSOCIATED(integrals%hdkh)) THEN
1188         DEALLOCATE (integrals%hdkh)
1189      END IF
1190
1191      integrals%zorastat = 0
1192      integrals%dkhstat = 0
1193
1194   END SUBROUTINE atom_relint_release
1195
1196! **************************************************************************************************
1197!> \brief Calculate model potential. It is useful to guess atomic orbitals.
1198!> \param modpot  model potential
1199!> \param grid    atomic radial grid
1200!> \param zcore   nuclear charge
1201! **************************************************************************************************
1202   SUBROUTINE calculate_model_potential(modpot, grid, zcore)
1203      REAL(dp), DIMENSION(:), INTENT(INOUT)              :: modpot
1204      TYPE(grid_atom_type), INTENT(IN)                   :: grid
1205      REAL(dp), INTENT(IN)                               :: zcore
1206
1207      CHARACTER(len=*), PARAMETER :: routineN = 'calculate_model_potential', &
1208         routineP = moduleN//':'//routineN
1209
1210      INTEGER                                            :: i, ii, l, ll, n, z
1211      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: pot, rho
1212      TYPE(atom_state)                                   :: state
1213
1214      n = SIZE(modpot)
1215      ALLOCATE (rho(n), pot(n))
1216
1217      ! fill default occupation
1218      state%core = 0._dp
1219      state%occ = 0._dp
1220      state%multiplicity = -1
1221      z = NINT(zcore)
1222      DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
1223         IF (ptable(z)%e_conv(l) /= 0) THEN
1224            state%maxl_occ = l
1225            ll = 2*(2*l + 1)
1226            DO i = 1, SIZE(state%occ, 2)
1227               ii = ptable(z)%e_conv(l) - (i - 1)*ll
1228               IF (ii <= ll) THEN
1229                  state%occ(l, i) = ii
1230                  EXIT
1231               ELSE
1232                  state%occ(l, i) = ll
1233               END IF
1234            END DO
1235         END IF
1236      END DO
1237
1238      modpot = -zcore/grid%rad(:)
1239
1240      ! Coulomb potential
1241      CALL slater_density(rho, pot, NINT(zcore), state, grid)
1242      CALL coulomb_potential_numeric(pot, rho, grid)
1243      modpot = modpot + pot
1244
1245      ! XC potential
1246      CALL wigner_slater_functional(rho, pot)
1247      modpot = modpot + pot
1248
1249      DEALLOCATE (rho, pot)
1250
1251   END SUBROUTINE calculate_model_potential
1252
1253END MODULE atom_operators
1254