LCOV - code coverage report
Current view: top level - src - auto_basis.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 68.1 % 351 239
Test Date: 2025-07-25 12:55:17 Functions: 62.5 % 8 5

            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   Automatic generation of auxiliary basis sets of different kind
      10              : !> \author  JGH
      11              : !>
      12              : !> <b>Modification history:</b>
      13              : !> - 11.2017 creation [JGH]
      14              : ! **************************************************************************************************
      15              : MODULE auto_basis
      16              :    USE aux_basis_set,                   ONLY: create_aux_basis
      17              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      18              :                                               gto_basis_set_type,&
      19              :                                               sort_gto_basis_set
      20              :    USE bibliography,                    ONLY: Stoychev2016,&
      21              :                                               cite_reference
      22              :    USE kinds,                           ONLY: default_string_length,&
      23              :                                               dp
      24              :    USE mathconstants,                   ONLY: dfac,&
      25              :                                               fac,&
      26              :                                               gamma1,&
      27              :                                               pi,&
      28              :                                               rootpi
      29              :    USE orbital_pointers,                ONLY: init_orbital_pointers
      30              :    USE periodic_table,                  ONLY: get_ptable_info
      31              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      32              :                                               qs_kind_type
      33              : #include "./base/base_uses.f90"
      34              : 
      35              :    IMPLICIT NONE
      36              : 
      37              :    PRIVATE
      38              : 
      39              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'auto_basis'
      40              : 
      41              :    PUBLIC :: create_ri_aux_basis_set, create_lri_aux_basis_set, &
      42              :              create_oce_basis
      43              : 
      44              : CONTAINS
      45              : 
      46              : ! **************************************************************************************************
      47              : !> \brief Create a RI_AUX basis set using some heuristics
      48              : !> \param ri_aux_basis_set ...
      49              : !> \param qs_kind ...
      50              : !> \param basis_cntrl ...
      51              : !> \param basis_type ...
      52              : !> \param basis_sort ...
      53              : !> \date    01.11.2017
      54              : !> \author  JGH
      55              : ! **************************************************************************************************
      56          292 :    SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis_type, basis_sort)
      57              :       TYPE(gto_basis_set_type), POINTER                  :: ri_aux_basis_set
      58              :       TYPE(qs_kind_type), INTENT(IN)                     :: qs_kind
      59              :       INTEGER, INTENT(IN)                                :: basis_cntrl
      60              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: basis_type
      61              :       INTEGER, INTENT(IN), OPTIONAL                      :: basis_sort
      62              : 
      63              :       CHARACTER(LEN=2)                                   :: element_symbol
      64              :       CHARACTER(LEN=default_string_length)               :: bsname
      65              :       INTEGER                                            :: i, j, jj, l, laux, linc, lmax, lval, lx, &
      66              :                                                             nsets, nx, z
      67              :       INTEGER, DIMENSION(0:18)                           :: nval
      68              :       INTEGER, DIMENSION(0:9, 1:20)                      :: nl
      69              :       INTEGER, DIMENSION(1:3)                            :: ls1, ls2, npgf
      70          292 :       INTEGER, DIMENSION(:), POINTER                     :: econf
      71              :       REAL(KIND=dp)                                      :: xv, zval
      72          292 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: zet
      73              :       REAL(KIND=dp), DIMENSION(0:18)                     :: bv, bval, fv, peff, pend, pmax, pmin
      74              :       REAL(KIND=dp), DIMENSION(0:9)                      :: zeff, zmax, zmin
      75              :       REAL(KIND=dp), DIMENSION(3)                        :: amax, amin, bmin
      76              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
      77              : 
      78              :       !
      79          292 :       CALL cite_reference(Stoychev2016)
      80              :       !
      81              :       bv(0:18) = (/1.8_dp, 2.0_dp, 2.2_dp, 2.2_dp, 2.3_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, &
      82          292 :                    3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp/)
      83              :       fv(0:18) = (/20.0_dp, 4.0_dp, 4.0_dp, 3.5_dp, 2.5_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, &
      84          292 :                    2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp/)
      85              :       !
      86          292 :       CPASSERT(.NOT. ASSOCIATED(ri_aux_basis_set))
      87          292 :       NULLIFY (orb_basis_set)
      88          292 :       IF (.NOT. PRESENT(basis_type)) THEN
      89          234 :          CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
      90              :       ELSE
      91           58 :          CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type=basis_type)
      92              :       END IF
      93          292 :       IF (ASSOCIATED(orb_basis_set)) THEN
      94              :          ! BASIS_SET ORB NONE associates the pointer orb_basis_set, but does not contain
      95              :          ! any actual basis functions. Therefore, we catch it here to avoid spurious autogenerated
      96              :          ! RI_AUX basis sets.
      97         1310 :          IF (SUM(orb_basis_set%nsgf_set) == 0) THEN
      98              :             CALL cp_abort(__LOCATION__, &
      99              :                           "Cannot autocreate RI_AUX basis set for at least one of the given "// &
     100              :                           "primary basis sets due to missing exponents. If you have invoked BASIS_SET NONE, "// &
     101            0 :                           "you should state BASIS_SET RI_AUX NONE explicitly in the input.")
     102              :          END IF
     103          292 :          CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
     104              :          !Note: RI basis coud require lmax up to 2*orb_lmax. This ensures that all orbital pointers
     105              :          !      are properly initialized before building the basis
     106          292 :          CALL init_orbital_pointers(2*lmax)
     107          292 :          CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
     108          292 :          CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
     109          292 :          CALL get_ptable_info(element_symbol, ielement=z)
     110          292 :          lval = 0
     111         1436 :          DO l = 0, MAXVAL(UBOUND(econf))
     112         1144 :             IF (econf(l) > 0) lval = l
     113              :          END DO
     114         1144 :          IF (SUM(econf) /= NINT(zval)) THEN
     115            0 :             CPWARN("Valence charge and electron configuration not consistent")
     116              :          END IF
     117          292 :          pend = 0.0_dp
     118          292 :          linc = 1
     119          292 :          IF (z > 18) linc = 2
     120          440 :          SELECT CASE (basis_cntrl)
     121              :          CASE (0)
     122          148 :             laux = MAX(2*lval, lmax + linc)
     123              :          CASE (1)
     124          140 :             laux = MAX(2*lval, lmax + linc)
     125              :          CASE (2)
     126            0 :             laux = MAX(2*lval, lmax + linc + 1)
     127              :          CASE (3)
     128            4 :             laux = MAX(2*lmax, lmax + linc + 2)
     129              :          CASE DEFAULT
     130          292 :             CPABORT("Invalid value of control variable")
     131              :          END SELECT
     132              :          !
     133          346 :          DO l = 2*lmax + 1, laux
     134           54 :             xv = peff(2*lmax)
     135           54 :             pmin(l) = xv
     136           54 :             pmax(l) = xv
     137           54 :             peff(l) = xv
     138          346 :             pend(l) = xv
     139              :          END DO
     140              :          !
     141         1256 :          DO l = 0, laux
     142          964 :             IF (l <= 2*lval) THEN
     143          660 :                pend(l) = MIN(fv(l)*peff(l), pmax(l))
     144          660 :                bval(l) = 1.8_dp
     145              :             ELSE
     146          304 :                pend(l) = peff(l)
     147          304 :                bval(l) = bv(l)
     148              :             END IF
     149          964 :             xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
     150         1256 :             nval(l) = MAX(CEILING(xv), 0)
     151              :          END DO
     152              :          ! first set include valence only
     153          292 :          nsets = 1
     154          292 :          ls1(1) = 0
     155          292 :          ls2(1) = lval
     156          442 :          DO l = lval + 1, laux
     157          416 :             IF (nval(l) < nval(lval) - 1) EXIT
     158          442 :             ls2(1) = l
     159              :          END DO
     160              :          ! second set up to 2*lval
     161          292 :          IF (laux > ls2(1)) THEN
     162          266 :             IF (lval == 0 .OR. 2*lval <= ls2(1) + 1) THEN
     163          266 :                nsets = 2
     164          266 :                ls1(2) = ls2(1) + 1
     165          266 :                ls2(2) = laux
     166              :             ELSE
     167            0 :                nsets = 2
     168            0 :                ls1(2) = ls2(1) + 1
     169            0 :                ls2(2) = MIN(2*lval, laux)
     170            0 :                lx = ls2(2)
     171            0 :                DO l = lx + 1, laux
     172            0 :                   IF (nval(l) < nval(lx) - 1) EXIT
     173            0 :                   ls2(2) = l
     174              :                END DO
     175            0 :                IF (laux > ls2(2)) THEN
     176            0 :                   nsets = 3
     177            0 :                   ls1(3) = ls2(2) + 1
     178            0 :                   ls2(3) = laux
     179              :                END IF
     180              :             END IF
     181              :          END IF
     182              :          !
     183          292 :          amax = 0.0
     184         1168 :          amin = HUGE(0.0_dp)
     185         1168 :          bmin = HUGE(0.0_dp)
     186          850 :          DO i = 1, nsets
     187         1522 :             DO j = ls1(i), ls2(i)
     188          964 :                amax(i) = MAX(amax(i), pend(j))
     189          964 :                amin(i) = MIN(amin(i), pmin(j))
     190         1522 :                bmin(i) = MIN(bmin(i), bval(j))
     191              :             END DO
     192          558 :             xv = LOG(amax(i)/amin(i))/LOG(bmin(i)) + 1.e-10_dp
     193          850 :             npgf(i) = MAX(CEILING(xv), 0)
     194              :          END DO
     195          850 :          nx = MAXVAL(npgf(1:nsets))
     196         1168 :          ALLOCATE (zet(nx, nsets))
     197         6018 :          zet = 0.0_dp
     198          292 :          nl = 0
     199          850 :          DO i = 1, nsets
     200         3712 :             DO j = 1, npgf(i)
     201         3154 :                jj = npgf(i) - j + 1
     202         3712 :                zet(jj, i) = amin(i)*bmin(i)**(j - 1)
     203              :             END DO
     204         1814 :             DO l = ls1(i), ls2(i)
     205         1522 :                nl(l, i) = nval(l)
     206              :             END DO
     207              :          END DO
     208          292 :          bsname = TRIM(element_symbol)//"-RI-AUX-"//TRIM(orb_basis_set%name)
     209              :          !
     210          292 :          CALL create_aux_basis(ri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
     211              : 
     212          292 :          DEALLOCATE (zet)
     213              : 
     214          876 :          IF (PRESENT(basis_sort)) THEN
     215          186 :             CALL sort_gto_basis_set(ri_aux_basis_set, basis_sort)
     216              :          END IF
     217              : 
     218              :       END IF
     219              : 
     220          584 :    END SUBROUTINE create_ri_aux_basis_set
     221              : ! **************************************************************************************************
     222              : !> \brief Create a LRI_AUX basis set using some heuristics
     223              : !> \param lri_aux_basis_set ...
     224              : !> \param qs_kind ...
     225              : !> \param basis_cntrl ...
     226              : !> \param exact_1c_terms ...
     227              : !> \param tda_kernel ...
     228              : !> \date    01.11.2017
     229              : !> \author  JGH
     230              : ! **************************************************************************************************
     231           30 :    SUBROUTINE create_lri_aux_basis_set(lri_aux_basis_set, qs_kind, basis_cntrl, &
     232              :                                        exact_1c_terms, tda_kernel)
     233              :       TYPE(gto_basis_set_type), POINTER                  :: lri_aux_basis_set
     234              :       TYPE(qs_kind_type), INTENT(IN)                     :: qs_kind
     235              :       INTEGER, INTENT(IN)                                :: basis_cntrl
     236              :       LOGICAL, INTENT(IN), OPTIONAL                      :: exact_1c_terms, tda_kernel
     237              : 
     238              :       CHARACTER(LEN=2)                                   :: element_symbol
     239              :       CHARACTER(LEN=default_string_length)               :: bsname
     240              :       INTEGER                                            :: i, j, l, laux, linc, lm, lmax, lval, n1, &
     241              :                                                             n2, nsets, z
     242              :       INTEGER, DIMENSION(0:18)                           :: nval
     243              :       INTEGER, DIMENSION(0:9, 1:50)                      :: nl
     244              :       INTEGER, DIMENSION(1:50)                           :: ls1, ls2, npgf
     245           30 :       INTEGER, DIMENSION(:), POINTER                     :: econf
     246              :       LOGICAL                                            :: e1terms, kernel_basis
     247              :       REAL(KIND=dp)                                      :: xv, zval
     248           30 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: zet
     249              :       REAL(KIND=dp), DIMENSION(0:18)                     :: bval, peff, pend, pmax, pmin
     250              :       REAL(KIND=dp), DIMENSION(0:9)                      :: zeff, zmax, zmin
     251              :       REAL(KIND=dp), DIMENSION(4)                        :: bv, bx
     252              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     253              : 
     254              :       !
     255           30 :       IF (PRESENT(exact_1c_terms)) THEN
     256           30 :          e1terms = exact_1c_terms
     257              :       ELSE
     258              :          e1terms = .FALSE.
     259              :       END IF
     260           30 :       IF (PRESENT(tda_kernel)) THEN
     261           12 :          kernel_basis = tda_kernel
     262              :       ELSE
     263              :          kernel_basis = .FALSE.
     264              :       END IF
     265           30 :       IF (kernel_basis .AND. e1terms) THEN
     266            0 :          CALL cp_warn(__LOCATION__, "LRI Kernel basis generation will ignore exact 1C term option.")
     267              :       END IF
     268              :       !
     269           30 :       CPASSERT(.NOT. ASSOCIATED(lri_aux_basis_set))
     270           30 :       NULLIFY (orb_basis_set)
     271           30 :       CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
     272           30 :       IF (ASSOCIATED(orb_basis_set)) THEN
     273           30 :          CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
     274           30 :          CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
     275           30 :          CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
     276           30 :          CALL get_ptable_info(element_symbol, ielement=z)
     277           30 :          lval = 0
     278          106 :          DO l = 0, MAXVAL(UBOUND(econf))
     279           76 :             IF (econf(l) > 0) lval = l
     280              :          END DO
     281           76 :          IF (SUM(econf) /= NINT(zval)) THEN
     282            0 :             CPWARN("Valence charge and electron configuration not consistent")
     283              :          END IF
     284              :          !
     285           30 :          linc = 1
     286           30 :          IF (z > 18) linc = 2
     287           30 :          pend = 0.0_dp
     288           30 :          IF (kernel_basis) THEN
     289           12 :             bv(1:4) = (/3.20_dp, 2.80_dp, 2.40_dp, 2.00_dp/)
     290           12 :             bx(1:4) = (/4.00_dp, 3.50_dp, 3.00_dp, 2.50_dp/)
     291              :             !
     292           12 :             SELECT CASE (basis_cntrl)
     293              :             CASE (0)
     294            0 :                laux = lval + 1
     295              :             CASE (1)
     296           12 :                laux = MAX(lval + 1, lmax)
     297              :             CASE (2)
     298            0 :                laux = MAX(lval + 2, lmax + 1)
     299              :             CASE (3)
     300            0 :                laux = MAX(lval + 3, lmax + 2)
     301            0 :                laux = MIN(laux, 2 + linc)
     302              :             CASE DEFAULT
     303           12 :                CPABORT("Invalid value of control variable")
     304              :             END SELECT
     305              :          ELSE
     306           18 :             bv(1:4) = (/2.00_dp, 1.90_dp, 1.80_dp, 1.80_dp/)
     307           18 :             bx(1:4) = (/2.60_dp, 2.40_dp, 2.20_dp, 2.20_dp/)
     308              :             !
     309           18 :             SELECT CASE (basis_cntrl)
     310              :             CASE (0)
     311            0 :                laux = MAX(2*lval, lmax + linc)
     312            0 :                laux = MIN(laux, 2 + linc)
     313              :             CASE (1)
     314           18 :                laux = MAX(2*lval, lmax + linc)
     315           18 :                laux = MIN(laux, 3 + linc)
     316              :             CASE (2)
     317            0 :                laux = MAX(2*lval, lmax + linc + 1)
     318            0 :                laux = MIN(laux, 4 + linc)
     319              :             CASE (3)
     320            0 :                laux = MAX(2*lval, lmax + linc + 1)
     321            0 :                laux = MIN(laux, 4 + linc)
     322              :             CASE DEFAULT
     323           18 :                CPABORT("Invalid value of control variable")
     324              :             END SELECT
     325              :          END IF
     326              :          !
     327           30 :          DO l = 2*lmax + 1, laux
     328            0 :             pmin(l) = pmin(2*lmax)
     329            0 :             pmax(l) = pmax(2*lmax)
     330           30 :             peff(l) = peff(2*lmax)
     331              :          END DO
     332              :          !
     333           30 :          nval = 0
     334           30 :          IF (exact_1c_terms) THEN
     335            0 :             DO l = 0, laux
     336            0 :                IF (l <= lval + 1) THEN
     337            0 :                   pend(l) = zmax(l) + 1.0_dp
     338            0 :                   bval(l) = bv(basis_cntrl + 1)
     339              :                ELSE
     340            0 :                   pend(l) = 2.0_dp*peff(l)
     341            0 :                   bval(l) = bx(basis_cntrl + 1)
     342              :                END IF
     343            0 :                pmin(l) = zmin(l)
     344            0 :                xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
     345            0 :                nval(l) = MAX(CEILING(xv), 0)
     346            0 :                bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
     347              :             END DO
     348              :          ELSE
     349          124 :             DO l = 0, laux
     350           94 :                IF (l <= lval + 1) THEN
     351           76 :                   pend(l) = pmax(l)
     352           76 :                   bval(l) = bv(basis_cntrl + 1)
     353           76 :                   pmin(l) = zmin(l)
     354              :                ELSE
     355           18 :                   pend(l) = 4.0_dp*peff(l)
     356           18 :                   bval(l) = bx(basis_cntrl + 1)
     357              :                END IF
     358           94 :                xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
     359           94 :                nval(l) = MAX(CEILING(xv), 0)
     360          124 :                bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
     361              :             END DO
     362              :          END IF
     363              :          !
     364           30 :          lm = MIN(2*lval, 3)
     365           92 :          n1 = MAXVAL(nval(0:lm))
     366           30 :          IF (laux < lm + 1) THEN
     367              :             n2 = 0
     368              :          ELSE
     369           56 :             n2 = MAXVAL(nval(lm + 1:laux))
     370              :          END IF
     371              :          !
     372           30 :          nsets = n1 + n2
     373           90 :          ALLOCATE (zet(1, nsets))
     374          830 :          zet = 0.0_dp
     375           30 :          nl = 0
     376          152 :          j = MAXVAL(MAXLOC(nval(0:lm)))
     377          280 :          DO i = 1, n1
     378          250 :             ls1(i) = 0
     379          250 :             ls2(i) = lm
     380          250 :             npgf(i) = 1
     381          250 :             zet(1, i) = pmin(j)*bval(j)**(i - 1)
     382          786 :             DO l = 0, lm
     383          756 :                nl(l, i) = 1
     384              :             END DO
     385              :          END DO
     386           30 :          j = lm + 1
     387          180 :          DO i = n1 + 1, nsets
     388          150 :             ls1(i) = lm + 1
     389          150 :             ls2(i) = laux
     390          150 :             npgf(i) = 1
     391          150 :             zet(1, i) = pmin(j)*bval(j)**(i - n1 - 1)
     392          418 :             DO l = lm + 1, laux
     393          388 :                nl(l, i) = 1
     394              :             END DO
     395              :          END DO
     396              :          !
     397           30 :          bsname = TRIM(element_symbol)//"-LRI-AUX-"//TRIM(orb_basis_set%name)
     398              :          !
     399           30 :          CALL create_aux_basis(lri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
     400              :          !
     401           90 :          DEALLOCATE (zet)
     402              :       END IF
     403              : 
     404           60 :    END SUBROUTINE create_lri_aux_basis_set
     405              : 
     406              : ! **************************************************************************************************
     407              : !> \brief ...
     408              : !> \param oce_basis ...
     409              : !> \param orb_basis ...
     410              : !> \param lmax_oce ...
     411              : !> \param nbas_oce ...
     412              : ! **************************************************************************************************
     413           12 :    SUBROUTINE create_oce_basis(oce_basis, orb_basis, lmax_oce, nbas_oce)
     414              :       TYPE(gto_basis_set_type), POINTER                  :: oce_basis, orb_basis
     415              :       INTEGER, INTENT(IN)                                :: lmax_oce, nbas_oce
     416              : 
     417              :       CHARACTER(LEN=default_string_length)               :: bsname
     418              :       INTEGER                                            :: i, l, lmax, lx, nset, nx
     419           12 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: lmin, lset, npgf
     420           12 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: nl
     421           12 :       INTEGER, DIMENSION(:), POINTER                     :: npgf_orb
     422              :       REAL(KIND=dp)                                      :: cval, x, z0, z1
     423           12 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: zet
     424              :       REAL(KIND=dp), DIMENSION(0:9)                      :: zeff, zmax, zmin
     425              : 
     426           12 :       CALL get_basis_keyfigures(orb_basis, lmax, zmin, zmax, zeff)
     427           12 :       IF (nbas_oce < 1) THEN
     428           12 :          CALL get_gto_basis_set(gto_basis_set=orb_basis, nset=nset, npgf=npgf_orb)
     429           54 :          nx = SUM(npgf_orb(1:nset))
     430              :       ELSE
     431              :          nx = 0
     432              :       END IF
     433           12 :       nset = MAX(nbas_oce, nx)
     434           12 :       lx = MAX(lmax_oce, lmax)
     435              :       !
     436           12 :       bsname = "OCE-"//TRIM(orb_basis%name)
     437          108 :       ALLOCATE (lmin(nset), lset(nset), nl(0:9, nset), npgf(nset), zet(1, nset))
     438          114 :       lmin = 0
     439          114 :       lset = 0
     440         1134 :       nl = 1
     441          114 :       npgf = 1
     442          216 :       zet = 0.0_dp
     443              :       !
     444           56 :       z0 = MINVAL(zmin(0:lmax))
     445           56 :       z1 = MAXVAL(zmax(0:lmax))
     446           12 :       x = 1.0_dp/REAL(nset - 1, KIND=dp)
     447           12 :       cval = (z1/z0)**x
     448           12 :       zet(1, nset) = z0
     449          102 :       DO i = nset - 1, 1, -1
     450          102 :          zet(1, i) = zet(1, i + 1)*cval
     451              :       END DO
     452          114 :       DO i = 1, nset
     453          102 :          x = zet(1, i)
     454          284 :          DO l = 1, lmax
     455          182 :             z1 = 1.05_dp*zmax(l)
     456          284 :             IF (x < z1) lset(i) = l
     457              :          END DO
     458          114 :          IF (lset(i) == lmax) lset(i) = lx
     459              :       END DO
     460              :       !
     461           12 :       CALL create_aux_basis(oce_basis, bsname, nset, lmin, lset, nl, npgf, zet)
     462              :       !
     463           12 :       DEALLOCATE (lmin, lset, nl, npgf, zet)
     464              : 
     465           12 :    END SUBROUTINE create_oce_basis
     466              : ! **************************************************************************************************
     467              : !> \brief ...
     468              : !> \param basis_set ...
     469              : !> \param lmax ...
     470              : !> \param zmin ...
     471              : !> \param zmax ...
     472              : !> \param zeff ...
     473              : ! **************************************************************************************************
     474          334 :    SUBROUTINE get_basis_keyfigures(basis_set, lmax, zmin, zmax, zeff)
     475              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
     476              :       INTEGER, INTENT(OUT)                               :: lmax
     477              :       REAL(KIND=dp), DIMENSION(0:9), INTENT(OUT)         :: zmin, zmax, zeff
     478              : 
     479              :       INTEGER                                            :: i, ipgf, iset, ishell, j, l, nset
     480          334 :       INTEGER, DIMENSION(:), POINTER                     :: lm, npgf, nshell
     481          334 :       INTEGER, DIMENSION(:, :), POINTER                  :: lshell
     482              :       REAL(KIND=dp)                                      :: aeff, gcca, gccb, kval, rexp, rint, rno, &
     483              :                                                             zeta
     484          334 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
     485          334 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
     486              : 
     487              :       CALL get_gto_basis_set(gto_basis_set=basis_set, &
     488              :                              nset=nset, &
     489              :                              nshell=nshell, &
     490              :                              npgf=npgf, &
     491              :                              l=lshell, &
     492              :                              lmax=lm, &
     493              :                              zet=zet, &
     494          334 :                              gcc=gcc)
     495              : 
     496         1432 :       lmax = MAXVAL(lm)
     497          334 :       CPASSERT(lmax <= 9)
     498              : 
     499          334 :       zmax = 0.0_dp
     500         3674 :       zmin = HUGE(0.0_dp)
     501          334 :       zeff = 0.0_dp
     502              : 
     503         1432 :       DO iset = 1, nset
     504              :          ! zmin zmax
     505         3564 :          DO ipgf = 1, npgf(iset)
     506         7602 :             DO ishell = 1, nshell(iset)
     507         4038 :                l = lshell(ishell, iset)
     508         4038 :                zeta = zet(ipgf, iset)
     509         4038 :                zmax(l) = MAX(zmax(l), zeta)
     510         6504 :                zmin(l) = MIN(zmin(l), zeta)
     511              :             END DO
     512              :          END DO
     513              :          ! zeff
     514         2932 :          DO ishell = 1, nshell(iset)
     515         1500 :             l = lshell(ishell, iset)
     516         1500 :             kval = fac(l + 1)**2*2._dp**(2*l + 1)/fac(2*l + 2)
     517         1500 :             rexp = 0.0_dp
     518         1500 :             rno = 0.0_dp
     519         5538 :             DO i = 1, npgf(iset)
     520         4038 :                gcca = gcc(i, ishell, iset)
     521        22036 :                DO j = 1, npgf(iset)
     522        16498 :                   zeta = zet(i, iset) + zet(j, iset)
     523        16498 :                   gccb = gcc(j, ishell, iset)
     524        16498 :                   rint = 0.5_dp*fac(l + 1)/zeta**(l + 2)
     525        16498 :                   rexp = rexp + gcca*gccb*rint
     526        16498 :                   rint = rootpi*0.5_dp**(l + 2)*dfac(2*l + 1)/zeta**(l + 1.5_dp)
     527        20536 :                   rno = rno + gcca*gccb*rint
     528              :                END DO
     529              :             END DO
     530         1500 :             rexp = rexp/rno
     531         1500 :             aeff = (fac(l + 1)/dfac(2*l + 1))**2*2._dp**(2*l + 1)/(pi*rexp**2)
     532         2598 :             zeff(l) = MAX(zeff(l), aeff)
     533              :          END DO
     534              :       END DO
     535              : 
     536          334 :    END SUBROUTINE get_basis_keyfigures
     537              : 
     538              : ! **************************************************************************************************
     539              : !> \brief ...
     540              : !> \param lmax ...
     541              : !> \param zmin ...
     542              : !> \param zmax ...
     543              : !> \param zeff ...
     544              : !> \param pmin ...
     545              : !> \param pmax ...
     546              : !> \param peff ...
     547              : ! **************************************************************************************************
     548          322 :    SUBROUTINE get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
     549              :       INTEGER, INTENT(IN)                                :: lmax
     550              :       REAL(KIND=dp), DIMENSION(0:9), INTENT(IN)          :: zmin, zmax, zeff
     551              :       REAL(KIND=dp), DIMENSION(0:18), INTENT(OUT)        :: pmin, pmax, peff
     552              : 
     553              :       INTEGER                                            :: l1, l2, la
     554              : 
     555         6440 :       pmin = HUGE(0.0_dp)
     556          322 :       pmax = 0.0_dp
     557          322 :       peff = 0.0_dp
     558              : 
     559         1062 :       DO l1 = 0, lmax
     560         2368 :          DO l2 = l1, lmax
     561         4788 :             DO la = l2 - l1, l2 + l1
     562         2742 :                pmax(la) = MAX(pmax(la), zmax(l1) + zmax(l2))
     563         2742 :                pmin(la) = MIN(pmin(la), zmin(l1) + zmin(l2))
     564         4048 :                peff(la) = MAX(peff(la), zeff(l1) + zeff(l2))
     565              :             END DO
     566              :          END DO
     567              :       END DO
     568              : 
     569          322 :    END SUBROUTINE get_basis_products
     570              : ! **************************************************************************************************
     571              : !> \brief ...
     572              : !> \param lm ...
     573              : !> \param npgf ...
     574              : !> \param nfun ...
     575              : !> \param zet ...
     576              : !> \param gcc ...
     577              : !> \param nfit ...
     578              : !> \param afit ...
     579              : !> \param amet ...
     580              : !> \param eval ...
     581              : ! **************************************************************************************************
     582            0 :    SUBROUTINE overlap_maximum(lm, npgf, nfun, zet, gcc, nfit, afit, amet, eval)
     583              :       INTEGER, INTENT(IN)                                :: lm, npgf, nfun
     584              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zet
     585              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: gcc
     586              :       INTEGER, INTENT(IN)                                :: nfit
     587              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: afit
     588              :       REAL(KIND=dp), INTENT(IN)                          :: amet
     589              :       REAL(KIND=dp), INTENT(OUT)                         :: eval
     590              : 
     591              :       INTEGER                                            :: i, ia, ib, info
     592              :       REAL(KIND=dp)                                      :: fij, fxij, intab, p, xij
     593            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: fx, tx, x2, xx
     594              : 
     595              :       ! SUM_i(fi M fi)
     596            0 :       fij = 0.0_dp
     597            0 :       DO ia = 1, npgf
     598            0 :          DO ib = 1, npgf
     599            0 :             p = zet(ia) + zet(ib) + amet
     600            0 :             intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
     601            0 :             DO i = 1, nfun
     602            0 :                fij = fij + gcc(ia, i)*gcc(ib, i)*intab
     603              :             END DO
     604              :          END DO
     605              :       END DO
     606              : 
     607              :       !Integrals (fi M xj)
     608            0 :       ALLOCATE (fx(nfit, nfun), tx(nfit, nfun))
     609            0 :       fx = 0.0_dp
     610            0 :       DO ia = 1, npgf
     611            0 :          DO ib = 1, nfit
     612            0 :             p = zet(ia) + afit(ib) + amet
     613            0 :             intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
     614            0 :             DO i = 1, nfun
     615            0 :                fx(ib, i) = fx(ib, i) + gcc(ia, i)*intab
     616              :             END DO
     617              :          END DO
     618              :       END DO
     619              : 
     620              :       !Integrals (xi M xj)
     621            0 :       ALLOCATE (xx(nfit, nfit), x2(nfit, nfit))
     622            0 :       DO ia = 1, nfit
     623            0 :          DO ib = 1, nfit
     624            0 :             p = afit(ia) + afit(ib) + amet
     625            0 :             xx(ia, ib) = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
     626              :          END DO
     627              :       END DO
     628              : 
     629              :       !Solve for tab
     630            0 :       tx(1:nfit, 1:nfun) = fx(1:nfit, 1:nfun)
     631            0 :       x2(1:nfit, 1:nfit) = xx(1:nfit, 1:nfit)
     632            0 :       CALL dposv("U", nfit, nfun, x2, nfit, tx, nfit, info)
     633            0 :       IF (info == 0) THEN
     634              :          ! value t*xx*t
     635              :          xij = 0.0_dp
     636            0 :          DO i = 1, nfun
     637            0 :             xij = xij + DOT_PRODUCT(tx(:, i), MATMUL(xx, tx(:, i)))
     638              :          END DO
     639              :          ! value t*fx
     640              :          fxij = 0.0_dp
     641            0 :          DO i = 1, nfun
     642            0 :             fxij = fxij + DOT_PRODUCT(tx(:, i), fx(:, i))
     643              :          END DO
     644              :          !
     645            0 :          eval = fij - 2.0_dp*fxij + xij
     646              :       ELSE
     647              :          ! error in solving for max overlap
     648            0 :          eval = 1.0e10_dp
     649              :       END IF
     650              : 
     651            0 :       DEALLOCATE (fx, xx, x2, tx)
     652              : 
     653            0 :    END SUBROUTINE overlap_maximum
     654              : ! **************************************************************************************************
     655              : !> \brief ...
     656              : !> \param x ...
     657              : !> \param n ...
     658              : !> \param eval ...
     659              : ! **************************************************************************************************
     660            0 :    SUBROUTINE neb_potential(x, n, eval)
     661              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: x
     662              :       INTEGER, INTENT(IN)                                :: n
     663              :       REAL(KIND=dp), INTENT(INOUT)                       :: eval
     664              : 
     665              :       INTEGER                                            :: i
     666              : 
     667            0 :       DO i = 2, n
     668            0 :          IF (x(i) < 1.5_dp) THEN
     669            0 :             eval = eval + 10.0_dp*(1.5_dp - x(i))**2
     670              :          END IF
     671              :       END DO
     672              : 
     673            0 :    END SUBROUTINE neb_potential
     674              : ! **************************************************************************************************
     675              : !> \brief ...
     676              : !> \param basis_set ...
     677              : !> \param lin ...
     678              : !> \param np ...
     679              : !> \param nf ...
     680              : !> \param zval ...
     681              : !> \param gcval ...
     682              : ! **************************************************************************************************
     683            0 :    SUBROUTINE get_basis_functions(basis_set, lin, np, nf, zval, gcval)
     684              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
     685              :       INTEGER, INTENT(IN)                                :: lin
     686              :       INTEGER, INTENT(OUT)                               :: np, nf
     687              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zval
     688              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gcval
     689              : 
     690              :       INTEGER                                            :: iset, ishell, j1, j2, jf, jp, l, nset
     691            0 :       INTEGER, DIMENSION(:), POINTER                     :: lm, npgf, nshell
     692            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: lshell
     693              :       LOGICAL                                            :: toadd
     694            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
     695            0 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
     696              : 
     697              :       CALL get_gto_basis_set(gto_basis_set=basis_set, &
     698              :                              nset=nset, &
     699              :                              nshell=nshell, &
     700              :                              npgf=npgf, &
     701              :                              l=lshell, &
     702              :                              lmax=lm, &
     703              :                              zet=zet, &
     704            0 :                              gcc=gcc)
     705              : 
     706            0 :       np = 0
     707            0 :       nf = 0
     708            0 :       DO iset = 1, nset
     709            0 :          toadd = .TRUE.
     710            0 :          DO ishell = 1, nshell(iset)
     711            0 :             l = lshell(ishell, iset)
     712            0 :             IF (l == lin) THEN
     713            0 :                nf = nf + 1
     714            0 :                IF (toadd) THEN
     715            0 :                   np = np + npgf(iset)
     716            0 :                   toadd = .FALSE.
     717              :                END IF
     718              :             END IF
     719              :          END DO
     720              :       END DO
     721            0 :       ALLOCATE (zval(np), gcval(np, nf))
     722            0 :       zval = 0.0_dp
     723            0 :       gcval = 0.0_dp
     724              :       !
     725              :       jp = 0
     726              :       jf = 0
     727            0 :       DO iset = 1, nset
     728            0 :          toadd = .TRUE.
     729            0 :          DO ishell = 1, nshell(iset)
     730            0 :             l = lshell(ishell, iset)
     731            0 :             IF (l == lin) THEN
     732            0 :                jf = jf + 1
     733            0 :                IF (toadd) THEN
     734            0 :                   j1 = jp + 1
     735            0 :                   j2 = jp + npgf(iset)
     736            0 :                   zval(j1:j2) = zet(1:npgf(iset), iset)
     737            0 :                   jp = jp + npgf(iset)
     738            0 :                   toadd = .FALSE.
     739              :                END IF
     740            0 :                gcval(j1:j2, jf) = gcc(1:npgf(iset), ishell, iset)
     741              :             END IF
     742              :          END DO
     743              :       END DO
     744              : 
     745            0 :    END SUBROUTINE get_basis_functions
     746              : 
     747            0 : END MODULE auto_basis
        

Generated by: LCOV version 2.0-1