LCOV - code coverage report
Current view: top level - src - lri_environment_init.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 77.2 % 302 233
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief initializes the environment for lri
      10              : !>        lri : local resolution of the identity
      11              : !> \par History
      12              : !>      created [06.2015]
      13              : !> \author Dorothea Golze
      14              : ! **************************************************************************************************
      15              : MODULE lri_environment_init
      16              :    USE ao_util,                         ONLY: exp_radius
      17              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18              :                                               get_atomic_kind_set
      19              :    USE basis_set_types,                 ONLY: copy_gto_basis_set,&
      20              :                                               gto_basis_set_type
      21              :    USE bibliography,                    ONLY: Golze2017a,&
      22              :                                               Golze2017b,&
      23              :                                               cite_reference
      24              :    USE cp_control_types,                ONLY: dft_control_type
      25              :    USE generic_shg_integrals,           ONLY: int_overlap_aba_shg
      26              :    USE generic_shg_integrals_init,      ONLY: contraction_matrix_shg,&
      27              :                                               contraction_matrix_shg_mix,&
      28              :                                               get_clebsch_gordon_coefficients
      29              :    USE input_section_types,             ONLY: section_vals_type,&
      30              :                                               section_vals_val_get
      31              :    USE kinds,                           ONLY: dp
      32              :    USE lri_environment_types,           ONLY: deallocate_bas_properties,&
      33              :                                               lri_env_create,&
      34              :                                               lri_environment_type
      35              :    USE mathconstants,                   ONLY: fac,&
      36              :                                               pi,&
      37              :                                               rootpi
      38              :    USE mathlib,                         ONLY: invert_matrix
      39              :    USE qs_environment_types,            ONLY: get_qs_env,&
      40              :                                               qs_environment_type
      41              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      42              :                                               qs_kind_type
      43              : #include "./base/base_uses.f90"
      44              : 
      45              :    IMPLICIT NONE
      46              : 
      47              :    PRIVATE
      48              : 
      49              : ! **************************************************************************************************
      50              : 
      51              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_environment_init'
      52              : 
      53              :    PUBLIC :: lri_env_init, lri_env_basis, lri_basis_init
      54              : 
      55              : ! **************************************************************************************************
      56              : 
      57              : CONTAINS
      58              : 
      59              : ! **************************************************************************************************
      60              : !> \brief initializes the lri env
      61              : !> \param lri_env ...
      62              : !> \param lri_section ...
      63              : ! **************************************************************************************************
      64          120 :    SUBROUTINE lri_env_init(lri_env, lri_section)
      65              : 
      66              :       TYPE(lri_environment_type), POINTER                :: lri_env
      67              :       TYPE(section_vals_type), POINTER                   :: lri_section
      68              : 
      69           60 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: radii
      70              : 
      71              :       NULLIFY (lri_env)
      72            0 :       ALLOCATE (lri_env)
      73           60 :       CALL lri_env_create(lri_env)
      74              : 
      75              :       ! init keywords
      76              :       ! use RI for local pp terms
      77              :       CALL section_vals_val_get(lri_section, "RI_STATISTIC", &
      78           60 :                                 l_val=lri_env%statistics)
      79              :       ! use exact one-center terms
      80              :       CALL section_vals_val_get(lri_section, "EXACT_1C_TERMS", &
      81           60 :                                 l_val=lri_env%exact_1c_terms)
      82              :       ! use RI for local pp terms
      83              :       CALL section_vals_val_get(lri_section, "PPL_RI", &
      84           60 :                                 l_val=lri_env%ppl_ri)
      85              :       ! check for debug (OS scheme)
      86              :       CALL section_vals_val_get(lri_section, "DEBUG_LRI_INTEGRALS", &
      87           60 :                                 l_val=lri_env%debug)
      88              :       ! integrals based on solid harmonic Gaussians
      89              :       CALL section_vals_val_get(lri_section, "SHG_LRI_INTEGRALS", &
      90           60 :                                 l_val=lri_env%use_shg_integrals)
      91              :       ! how to calculate inverse/pseuodinverse of overlap
      92              :       CALL section_vals_val_get(lri_section, "LRI_OVERLAP_MATRIX", &
      93           60 :                                 i_val=lri_env%lri_overlap_inv)
      94              :       CALL section_vals_val_get(lri_section, "MAX_CONDITION_NUM", &
      95           60 :                                 r_val=lri_env%cond_max)
      96              :       ! integrals threshold (aba, abb)
      97              :       CALL section_vals_val_get(lri_section, "EPS_O3_INT", &
      98           60 :                                 r_val=lri_env%eps_o3_int)
      99              :       ! RI SINV
     100              :       CALL section_vals_val_get(lri_section, "RI_SINV", &
     101           60 :                                 c_val=lri_env%ri_sinv_app)
     102              :       ! Distant Pair Approximation
     103              :       CALL section_vals_val_get(lri_section, "DISTANT_PAIR_APPROXIMATION", &
     104           60 :                                 l_val=lri_env%distant_pair_approximation)
     105              :       CALL section_vals_val_get(lri_section, "DISTANT_PAIR_METHOD", &
     106           60 :                                 c_val=lri_env%distant_pair_method)
     107           60 :       CALL section_vals_val_get(lri_section, "DISTANT_PAIR_RADII", r_vals=radii)
     108           60 :       CPASSERT(SIZE(radii) == 2)
     109           60 :       CPASSERT(radii(2) > radii(1))
     110           60 :       CPASSERT(radii(1) > 0.0_dp)
     111           60 :       lri_env%r_in = radii(1)
     112           60 :       lri_env%r_out = radii(2)
     113              : 
     114           60 :       CALL cite_reference(Golze2017b)
     115           60 :       IF (lri_env%use_shg_integrals) CALL cite_reference(Golze2017a)
     116              : 
     117           60 :    END SUBROUTINE lri_env_init
     118              : ! **************************************************************************************************
     119              : !> \brief initializes the lri env
     120              : !> \param ri_type ...
     121              : !> \param qs_env ...
     122              : !> \param lri_env ...
     123              : !> \param qs_kind_set ...
     124              : ! **************************************************************************************************
     125           60 :    SUBROUTINE lri_env_basis(ri_type, qs_env, lri_env, qs_kind_set)
     126              : 
     127              :       CHARACTER(len=*), INTENT(IN)                       :: ri_type
     128              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     129              :       TYPE(lri_environment_type), POINTER                :: lri_env
     130              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     131              : 
     132              :       INTEGER :: i, i1, i2, iat, ikind, ip, ipgf, iset, ishell, jp, l, lmax_ikind_orb, &
     133              :          lmax_ikind_ri, maxl_orb, maxl_ri, n1, n2, natom, nbas, nkind, nribas, nspin
     134           60 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     135              :       REAL(KIND=dp)                                      :: gcc, rad, rai, raj, xradius, zeta
     136           60 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: int_aux, norm
     137           60 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     138              :       TYPE(dft_control_type), POINTER                    :: dft_control
     139              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set, ri_basis_set
     140              : 
     141              :       ! initialize the basic basis sets (orb and ri)
     142           60 :       CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
     143           60 :       nkind = SIZE(atomic_kind_set)
     144          460 :       ALLOCATE (lri_env%orb_basis(nkind), lri_env%ri_basis(nkind))
     145           60 :       maxl_orb = 0
     146           60 :       maxl_ri = 0
     147          170 :       DO ikind = 1, nkind
     148          110 :          NULLIFY (orb_basis_set, ri_basis_set)
     149          110 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
     150          110 :          IF (ri_type == "LRI") THEN
     151           90 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="LRI_AUX")
     152           20 :          ELSE IF (ri_type == "P_LRI") THEN
     153           20 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="P_LRI_AUX")
     154            0 :          ELSE IF (ri_type == "RI") THEN
     155            0 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC")
     156              :          ELSE
     157            0 :             CPABORT('ri_type')
     158              :          END IF
     159          110 :          NULLIFY (lri_env%orb_basis(ikind)%gto_basis_set)
     160          110 :          NULLIFY (lri_env%ri_basis(ikind)%gto_basis_set)
     161          110 :          IF (ASSOCIATED(orb_basis_set)) THEN
     162          110 :             CALL copy_gto_basis_set(orb_basis_set, lri_env%orb_basis(ikind)%gto_basis_set)
     163          110 :             CALL copy_gto_basis_set(ri_basis_set, lri_env%ri_basis(ikind)%gto_basis_set)
     164              :          END IF
     165          264 :          lmax_ikind_orb = MAXVAL(lri_env%orb_basis(ikind)%gto_basis_set%lmax)
     166         1368 :          lmax_ikind_ri = MAXVAL(lri_env%ri_basis(ikind)%gto_basis_set%lmax)
     167          110 :          maxl_orb = MAX(maxl_orb, lmax_ikind_orb)
     168          170 :          maxl_ri = MAX(maxl_ri, lmax_ikind_ri)
     169              :       END DO
     170              : 
     171           60 :       IF ((ri_type == "LRI") .OR. (ri_type == "P_LRI")) THEN
     172              :          ! CG coefficients needed for lri integrals
     173           60 :          IF (ASSOCIATED(lri_env%cg_shg)) THEN
     174              :             CALL get_clebsch_gordon_coefficients(lri_env%cg_shg%cg_coeff, &
     175              :                                                  lri_env%cg_shg%cg_none0_list, &
     176              :                                                  lri_env%cg_shg%ncg_none0, &
     177           60 :                                                  maxl_orb, maxl_ri)
     178              :          END IF
     179           60 :          CALL lri_basis_init(lri_env)
     180              :          ! distant pair approximation
     181           60 :          IF (lri_env%distant_pair_approximation) THEN
     182              :             !
     183            0 :             SELECT CASE (lri_env%distant_pair_method)
     184              :             CASE ("EW")
     185              :                ! equal weight of 0.5
     186              :             CASE ("AW")
     187            0 :                ALLOCATE (lri_env%aradius(nkind))
     188            0 :                DO i = 1, nkind
     189            0 :                   orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
     190            0 :                   lri_env%aradius(i) = orb_basis_set%kind_radius
     191              :                END DO
     192              :             CASE ("SW")
     193            0 :                ALLOCATE (lri_env%wbas(nkind))
     194            0 :                DO i = 1, nkind
     195            0 :                   orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
     196            0 :                   n1 = orb_basis_set%nsgf
     197            0 :                   ALLOCATE (lri_env%wbas(i)%vec(n1))
     198            0 :                   DO iset = 1, orb_basis_set%nset
     199            0 :                      i1 = orb_basis_set%first_sgf(1, iset)
     200            0 :                      n2 = orb_basis_set%nshell(iset)
     201            0 :                      i2 = orb_basis_set%last_sgf(n2, iset)
     202            0 :                      lri_env%wbas(i)%vec(i1:i2) = orb_basis_set%set_radius(iset)
     203              :                   END DO
     204              :                END DO
     205              :             CASE ("LW")
     206           78 :                ALLOCATE (lri_env%wbas(nkind))
     207           46 :                DO i = 1, nkind
     208           30 :                   orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
     209           30 :                   n1 = orb_basis_set%nsgf
     210           90 :                   ALLOCATE (lri_env%wbas(i)%vec(n1))
     211           76 :                   DO iset = 1, orb_basis_set%nset
     212          182 :                      DO ishell = 1, orb_basis_set%nshell(iset)
     213          122 :                         i1 = orb_basis_set%first_sgf(ishell, iset)
     214          122 :                         i2 = orb_basis_set%last_sgf(ishell, iset)
     215          122 :                         l = orb_basis_set%l(ishell, iset)
     216          122 :                         xradius = 0.0_dp
     217          956 :                         DO ipgf = 1, orb_basis_set%npgf(iset)
     218          834 :                            gcc = orb_basis_set%gcc(ipgf, ishell, iset)
     219          834 :                            zeta = orb_basis_set%zet(ipgf, iset)
     220          834 :                            rad = exp_radius(l, zeta, 1.e-5_dp, gcc, rlow=xradius)
     221          956 :                            xradius = MAX(xradius, rad)
     222              :                         END DO
     223          430 :                         lri_env%wbas(i)%vec(i1:i2) = xradius
     224              :                      END DO
     225              :                   END DO
     226              :                END DO
     227              :             CASE DEFAULT
     228           18 :                CPABORT("Unknown DISTANT_PAIR_METHOD in LRI")
     229              :             END SELECT
     230              :             !
     231          172 :             ALLOCATE (lri_env%wmat(nkind, nkind))
     232           18 :             SELECT CASE (lri_env%distant_pair_method)
     233              :             CASE ("EW")
     234              :                ! equal weight of 0.5
     235            6 :                DO i1 = 1, nkind
     236            4 :                   n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
     237           14 :                   DO i2 = 1, nkind
     238            8 :                      n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
     239           32 :                      ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
     240          732 :                      lri_env%wmat(i1, i2)%mat(:, :) = 0.5_dp
     241              :                   END DO
     242              :                END DO
     243              :             CASE ("AW")
     244            0 :                DO i1 = 1, nkind
     245            0 :                   n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
     246            0 :                   DO i2 = 1, nkind
     247            0 :                      n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
     248            0 :                      ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
     249            0 :                      rai = lri_env%aradius(i1)**2
     250            0 :                      raj = lri_env%aradius(i2)**2
     251            0 :                      IF (raj > rai) THEN
     252            0 :                         lri_env%wmat(i1, i2)%mat(:, :) = 1.0_dp
     253              :                      ELSE
     254            0 :                         lri_env%wmat(i1, i2)%mat(:, :) = 0.0_dp
     255              :                      END IF
     256              :                   END DO
     257              :                END DO
     258              :             CASE ("SW", "LW")
     259           48 :                DO i1 = 1, nkind
     260           30 :                   n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
     261          104 :                   DO i2 = 1, nkind
     262           58 :                      n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
     263          232 :                      ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
     264          618 :                      DO ip = 1, SIZE(lri_env%wbas(i1)%vec)
     265          530 :                         rai = lri_env%wbas(i1)%vec(ip)**2
     266         5462 :                         DO jp = 1, SIZE(lri_env%wbas(i2)%vec)
     267         4874 :                            raj = lri_env%wbas(i2)%vec(jp)**2
     268         5404 :                            IF (raj > rai) THEN
     269         2000 :                               lri_env%wmat(i1, i2)%mat(ip, jp) = 1.0_dp
     270              :                            ELSE
     271         2874 :                               lri_env%wmat(i1, i2)%mat(ip, jp) = 0.0_dp
     272              :                            END IF
     273              :                         END DO
     274              :                      END DO
     275              :                   END DO
     276              :                END DO
     277              :             END SELECT
     278              :          END IF
     279            0 :       ELSE IF (ri_type == "RI") THEN
     280            0 :          ALLOCATE (lri_env%ri_fit)
     281            0 :          NULLIFY (lri_env%ri_fit%nvec)
     282            0 :          NULLIFY (lri_env%ri_fit%bas_ptr)
     283            0 :          CALL get_qs_env(qs_env=qs_env, natom=natom)
     284              :          ! initialize pointers to RI basis vector
     285            0 :          ALLOCATE (lri_env%ri_fit%bas_ptr(2, natom))
     286            0 :          ALLOCATE (kind_of(natom))
     287            0 :          CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     288            0 :          nbas = 0
     289            0 :          DO iat = 1, natom
     290            0 :             ikind = kind_of(iat)
     291            0 :             nribas = lri_env%ri_basis(ikind)%gto_basis_set%nsgf
     292            0 :             lri_env%ri_fit%bas_ptr(1, iat) = nbas + 1
     293            0 :             lri_env%ri_fit%bas_ptr(2, iat) = nbas + nribas
     294            0 :             nbas = nbas + nribas
     295              :          END DO
     296              :          ! initialize vector t
     297            0 :          CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     298            0 :          nspin = dft_control%nspins
     299            0 :          ALLOCATE (lri_env%ri_fit%tvec(nbas, nspin), lri_env%ri_fit%rm1t(nbas, nspin))
     300              :          ! initialize vector a, expansion of density
     301            0 :          ALLOCATE (lri_env%ri_fit%avec(nbas, nspin))
     302              :          ! initialize vector fout, R^(-1)*(f-p*n)
     303            0 :          ALLOCATE (lri_env%ri_fit%fout(nbas, nspin))
     304              :          ! initialize vector with RI basis integrated
     305            0 :          NULLIFY (norm, int_aux)
     306            0 :          nbas = lri_env%ri_fit%bas_ptr(2, natom)
     307            0 :          ALLOCATE (lri_env%ri_fit%nvec(nbas), lri_env%ri_fit%rm1n(nbas))
     308            0 :          ikind = 0
     309            0 :          DO iat = 1, natom
     310            0 :             IF (ikind /= kind_of(iat)) THEN
     311            0 :                ikind = kind_of(iat)
     312            0 :                ri_basis_set => lri_env%ri_basis(ikind)%gto_basis_set
     313            0 :                IF (ASSOCIATED(norm)) DEALLOCATE (norm)
     314            0 :                IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
     315            0 :                CALL basis_norm_s_func(ri_basis_set, norm)
     316            0 :                CALL basis_int(ri_basis_set, int_aux, norm)
     317              :             END IF
     318            0 :             nbas = SIZE(int_aux)
     319            0 :             i1 = lri_env%ri_fit%bas_ptr(1, iat)
     320            0 :             i2 = lri_env%ri_fit%bas_ptr(2, iat)
     321            0 :             lri_env%ri_fit%nvec(i1:i2) = int_aux(1:nbas)
     322              :          END DO
     323            0 :          IF (ASSOCIATED(norm)) DEALLOCATE (norm)
     324            0 :          IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
     325            0 :          DEALLOCATE (kind_of)
     326              :       ELSE
     327            0 :          CPABORT('ri_type')
     328              :       END IF
     329              : 
     330          120 :    END SUBROUTINE lri_env_basis
     331              : 
     332              : ! **************************************************************************************************
     333              : !> \brief initializes the lri basis: calculates the norm, self-overlap
     334              : !>        and integral of the ri basis
     335              : !> \param lri_env ...
     336              : ! **************************************************************************************************
     337           78 :    SUBROUTINE lri_basis_init(lri_env)
     338              :       TYPE(lri_environment_type), POINTER                :: lri_env
     339              : 
     340              :       INTEGER                                            :: ikind, nkind
     341           78 :       INTEGER, DIMENSION(:, :, :), POINTER               :: orb_index, ri_index
     342              :       REAL(KIND=dp)                                      :: delta
     343           78 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dovlp3
     344           78 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: orb_norm_r, ri_int_fbas, ri_norm_r, &
     345           78 :                                                             ri_norm_s
     346           78 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: orb_ovlp, ri_ovlp, ri_ovlp_inv
     347           78 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scon_orb, scon_ri
     348           78 :       REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: scon_mix
     349              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis, ri_basis
     350              : 
     351           78 :       IF (ASSOCIATED(lri_env)) THEN
     352           78 :          IF (ASSOCIATED(lri_env%orb_basis)) THEN
     353           78 :             CPASSERT(ASSOCIATED(lri_env%ri_basis))
     354           78 :             nkind = SIZE(lri_env%orb_basis)
     355           78 :             CALL deallocate_bas_properties(lri_env)
     356          362 :             ALLOCATE (lri_env%bas_prop(nkind))
     357          206 :             DO ikind = 1, nkind
     358          128 :                NULLIFY (orb_basis, ri_basis)
     359          128 :                orb_basis => lri_env%orb_basis(ikind)%gto_basis_set
     360          206 :                IF (ASSOCIATED(orb_basis)) THEN
     361          128 :                   ri_basis => lri_env%ri_basis(ikind)%gto_basis_set
     362          128 :                   CPASSERT(ASSOCIATED(ri_basis))
     363          128 :                   NULLIFY (ri_norm_r)
     364          128 :                   CALL basis_norm_radial(ri_basis, ri_norm_r)
     365          128 :                   NULLIFY (orb_norm_r)
     366          128 :                   CALL basis_norm_radial(orb_basis, orb_norm_r)
     367          128 :                   NULLIFY (ri_norm_s)
     368          128 :                   CALL basis_norm_s_func(ri_basis, ri_norm_s)
     369          128 :                   NULLIFY (ri_int_fbas)
     370          128 :                   CALL basis_int(ri_basis, ri_int_fbas, ri_norm_s)
     371          128 :                   lri_env%bas_prop(ikind)%int_fbas => ri_int_fbas
     372          128 :                   NULLIFY (ri_ovlp)
     373          128 :                   CALL basis_ovlp(ri_basis, ri_ovlp, ri_norm_r)
     374          128 :                   lri_env%bas_prop(ikind)%ri_ovlp => ri_ovlp
     375          128 :                   NULLIFY (orb_ovlp)
     376          128 :                   CALL basis_ovlp(orb_basis, orb_ovlp, orb_norm_r)
     377          128 :                   lri_env%bas_prop(ikind)%orb_ovlp => orb_ovlp
     378          128 :                   NULLIFY (scon_ri)
     379          128 :                   CALL contraction_matrix_shg(ri_basis, scon_ri)
     380          128 :                   lri_env%bas_prop(ikind)%scon_ri => scon_ri
     381          128 :                   NULLIFY (scon_orb)
     382          128 :                   CALL contraction_matrix_shg(orb_basis, scon_orb)
     383          128 :                   lri_env%bas_prop(ikind)%scon_orb => scon_orb
     384          128 :                   NULLIFY (scon_mix)
     385              :                   CALL contraction_matrix_shg_mix(orb_basis, ri_basis, &
     386          128 :                                                   orb_index, ri_index, scon_mix)
     387          128 :                   lri_env%bas_prop(ikind)%scon_mix => scon_mix
     388          128 :                   lri_env%bas_prop(ikind)%orb_index => orb_index
     389          128 :                   lri_env%bas_prop(ikind)%ri_index => ri_index
     390          640 :                   ALLOCATE (lri_env%bas_prop(ikind)%ovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf))
     391          768 :                   ALLOCATE (dovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf, 3))
     392              :                   CALL int_overlap_aba_shg(lri_env%bas_prop(ikind)%ovlp3, dovlp3, (/0.0_dp, 0.0_dp, 0.0_dp/), &
     393              :                                            orb_basis, orb_basis, ri_basis, scon_orb, &
     394              :                                            scon_mix, orb_index, ri_index, &
     395              :                                            lri_env%cg_shg%cg_coeff, &
     396              :                                            lri_env%cg_shg%cg_none0_list, &
     397              :                                            lri_env%cg_shg%ncg_none0, &
     398          128 :                                            calculate_forces=.FALSE.)
     399          128 :                   DEALLOCATE (orb_norm_r, ri_norm_r, ri_norm_s)
     400          128 :                   DEALLOCATE (dovlp3)
     401          512 :                   ALLOCATE (ri_ovlp_inv(ri_basis%nsgf, ri_basis%nsgf))
     402          128 :                   CALL invert_matrix(ri_ovlp, ri_ovlp_inv, delta, improve=.TRUE.)
     403          128 :                   lri_env%bas_prop(ikind)%ri_ovlp_inv => ri_ovlp_inv
     404              :                END IF
     405              :             END DO
     406              :          END IF
     407              :       END IF
     408              : 
     409          156 :    END SUBROUTINE lri_basis_init
     410              : 
     411              : ! **************************************************************************************************
     412              : !> \brief normalization for a contracted Gaussian s-function,
     413              : !>        spherical = cartesian Gaussian for s-functions
     414              : !> \param basis ...
     415              : !> \param norm ...
     416              : ! **************************************************************************************************
     417          128 :    SUBROUTINE basis_norm_s_func(basis, norm)
     418              : 
     419              :       TYPE(gto_basis_set_type), POINTER                  :: basis
     420              :       REAL(dp), DIMENSION(:), POINTER                    :: norm
     421              : 
     422              :       INTEGER                                            :: ipgf, iset, isgf, ishell, jpgf, l, nbas
     423              :       REAL(KIND=dp)                                      :: aai, aaj, cci, ccj, expa, ppl
     424              : 
     425          128 :       NULLIFY (norm)
     426              : 
     427          128 :       nbas = basis%nsgf
     428          384 :       ALLOCATE (norm(nbas))
     429        13778 :       norm = 0._dp
     430              : 
     431         1504 :       DO iset = 1, basis%nset
     432         5138 :          DO ishell = 1, basis%nshell(iset)
     433         3634 :             l = basis%l(ishell, iset)
     434         3634 :             IF (l /= 0) CYCLE
     435         1106 :             expa = 0.5_dp*REAL(2*l + 3, dp)
     436         1106 :             ppl = pi**(3._dp/2._dp)
     437         3588 :             DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
     438         2580 :                DO ipgf = 1, basis%npgf(iset)
     439         1474 :                   cci = basis%gcc(ipgf, ishell, iset)
     440         1474 :                   aai = basis%zet(ipgf, iset)
     441         6442 :                   DO jpgf = 1, basis%npgf(iset)
     442         3862 :                      ccj = basis%gcc(jpgf, ishell, iset)
     443         3862 :                      aaj = basis%zet(jpgf, iset)
     444         5336 :                      norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
     445              :                   END DO
     446              :                END DO
     447         4740 :                norm(isgf) = 1.0_dp/SQRT(norm(isgf))
     448              :             END DO
     449              :          END DO
     450              :       END DO
     451              : 
     452          128 :    END SUBROUTINE basis_norm_s_func
     453              : 
     454              : ! **************************************************************************************************
     455              : !> \brief normalization for radial part of contracted spherical Gaussian
     456              : !>        functions
     457              : !> \param basis ...
     458              : !> \param norm ...
     459              : ! **************************************************************************************************
     460          256 :    SUBROUTINE basis_norm_radial(basis, norm)
     461              : 
     462              :       TYPE(gto_basis_set_type), POINTER                  :: basis
     463              :       REAL(dp), DIMENSION(:), POINTER                    :: norm
     464              : 
     465              :       INTEGER                                            :: ipgf, iset, isgf, ishell, jpgf, l, nbas
     466              :       REAL(KIND=dp)                                      :: aai, aaj, cci, ccj, expa, ppl
     467              : 
     468          256 :       NULLIFY (norm)
     469              : 
     470          256 :       nbas = basis%nsgf
     471          768 :       ALLOCATE (norm(nbas))
     472        14884 :       norm = 0._dp
     473              : 
     474         1804 :       DO iset = 1, basis%nset
     475         5868 :          DO ishell = 1, basis%nshell(iset)
     476         4064 :             l = basis%l(ishell, iset)
     477         4064 :             expa = 0.5_dp*REAL(2*l + 3, dp)
     478         4064 :             ppl = fac(2*l + 2)*rootpi/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
     479        20240 :             DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
     480        36428 :                DO ipgf = 1, basis%npgf(iset)
     481        21800 :                   cci = basis%gcc(ipgf, ishell, iset)
     482        21800 :                   aai = basis%zet(ipgf, iset)
     483       105212 :                   DO jpgf = 1, basis%npgf(iset)
     484        68784 :                      ccj = basis%gcc(jpgf, ishell, iset)
     485        68784 :                      aaj = basis%zet(jpgf, iset)
     486        90584 :                      norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
     487              :                   END DO
     488              :                END DO
     489        18692 :                norm(isgf) = 1.0_dp/SQRT(norm(isgf))
     490              :             END DO
     491              :          END DO
     492              :       END DO
     493              : 
     494          256 :    END SUBROUTINE basis_norm_radial
     495              : 
     496              : !*****************************************************************************
     497              : !> \brief integral over a single (contracted) lri auxiliary basis function,
     498              : !>        integral is zero for all but s-functions
     499              : !> \param basis ...
     500              : !> \param int_aux ...
     501              : !> \param norm ...
     502              : ! **************************************************************************************************
     503          128 :    SUBROUTINE basis_int(basis, int_aux, norm)
     504              : 
     505              :       TYPE(gto_basis_set_type), POINTER                  :: basis
     506              :       REAL(dp), DIMENSION(:), POINTER                    :: int_aux, norm
     507              : 
     508              :       INTEGER                                            :: ipgf, iset, isgf, ishell, l, nbas
     509              :       REAL(KIND=dp)                                      :: aa, cc, pp
     510              : 
     511          128 :       nbas = basis%nsgf
     512          384 :       ALLOCATE (int_aux(nbas))
     513        13778 :       int_aux = 0._dp
     514              : 
     515         1504 :       DO iset = 1, basis%nset
     516         5138 :          DO ishell = 1, basis%nshell(iset)
     517         3634 :             l = basis%l(ishell, iset)
     518         3634 :             IF (l /= 0) CYCLE
     519         3588 :             DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
     520         6214 :                DO ipgf = 1, basis%npgf(iset)
     521         1474 :                   cc = basis%gcc(ipgf, ishell, iset)
     522         1474 :                   aa = basis%zet(ipgf, iset)
     523         1474 :                   pp = (pi/aa)**(3._dp/2._dp)
     524         2580 :                   int_aux(isgf) = int_aux(isgf) + norm(isgf)*cc*pp
     525              :                END DO
     526              :             END DO
     527              :          END DO
     528              :       END DO
     529              : 
     530          128 :    END SUBROUTINE basis_int
     531              : 
     532              : !*****************************************************************************
     533              : !> \brief self-overlap of lri basis for contracted spherical Gaussians.
     534              : !>        Overlap of radial part. Norm contains only normalization of radial
     535              : !>        part. Norm and overlap of spherical harmonics not explicitly
     536              : !>        calculated since this cancels for the self-overlap anyway.
     537              : !> \param basis ...
     538              : !> \param ovlp ...
     539              : !> \param norm ...
     540              : ! **************************************************************************************************
     541          256 :    SUBROUTINE basis_ovlp(basis, ovlp, norm)
     542              : 
     543              :       TYPE(gto_basis_set_type), POINTER                  :: basis
     544              :       REAL(dp), DIMENSION(:, :), POINTER                 :: ovlp
     545              :       REAL(dp), DIMENSION(:), POINTER                    :: norm
     546              : 
     547              :       INTEGER                                            :: ipgf, iset, isgf, ishell, jpgf, jset, &
     548              :                                                             jsgf, jshell, l, li, lj, m_i, m_j, nbas
     549              :       REAL(KIND=dp)                                      :: aai, aaj, cci, ccj, expa, norm_i, &
     550              :                                                             norm_j, oo, ppl
     551              : 
     552          256 :       nbas = basis%nsgf
     553         1024 :       ALLOCATE (ovlp(nbas, nbas))
     554      2343004 :       ovlp = 0._dp
     555              : 
     556         1804 :       DO iset = 1, basis%nset
     557         5868 :          DO ishell = 1, basis%nshell(iset)
     558         4064 :             li = basis%l(ishell, iset)
     559        52122 :             DO jset = 1, basis%nset
     560       192430 :                DO jshell = 1, basis%nshell(jset)
     561       141856 :                   lj = basis%l(jshell, jset)
     562       188366 :                   IF (li == lj) THEN
     563        34508 :                      l = li
     564        34508 :                      expa = 0.5_dp*REAL(2*l + 3, dp)
     565        34508 :                      ppl = fac(2*l + 2)*rootpi/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
     566       155708 :                      DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
     567       121200 :                         m_i = basis%m(isgf)
     568       775640 :                         DO jsgf = basis%first_sgf(jshell, jset), basis%last_sgf(jshell, jset)
     569       619932 :                            m_j = basis%m(jsgf)
     570       741132 :                            IF (m_i == m_j) THEN
     571       254220 :                               DO ipgf = 1, basis%npgf(iset)
     572       133020 :                                  cci = basis%gcc(ipgf, ishell, iset)
     573       133020 :                                  aai = basis%zet(ipgf, iset)
     574       133020 :                                  norm_i = norm(isgf)
     575       463872 :                                  DO jpgf = 1, basis%npgf(jset)
     576       209652 :                                     ccj = basis%gcc(jpgf, jshell, jset)
     577       209652 :                                     aaj = basis%zet(jpgf, jset)
     578       209652 :                                     oo = 1._dp/(aai + aaj)**expa
     579       209652 :                                     norm_j = norm(jsgf)
     580       342672 :                                     ovlp(isgf, jsgf) = ovlp(isgf, jsgf) + norm_i*norm_j*ppl*cci*ccj*oo
     581              :                                  END DO
     582              :                               END DO
     583              :                            END IF
     584              :                         END DO
     585              :                      END DO
     586              :                   END IF
     587              :                END DO
     588              :             END DO
     589              :          END DO
     590              :       END DO
     591              : 
     592          256 :    END SUBROUTINE basis_ovlp
     593              : 
     594              : END MODULE lri_environment_init
        

Generated by: LCOV version 2.0-1