LCOV - code coverage report
Current view: top level - src - atom_operators.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:f7c1873) Lines: 602 668 90.1 %
Date: 2024-03-26 06:29:31 Functions: 11 11 100.0 %

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

Generated by: LCOV version 1.15