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             : MODULE atom_set_basis
      10             :    USE ai_onecenter,                    ONLY: sg_overlap
      11             :    USE atom_types,                      ONLY: CGTO_BASIS,&
      12             :                                               GTO_BASIS,&
      13             :                                               atom_basis_type,&
      14             :                                               lmat
      15             :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      16             :                                               gto_basis_set_type
      17             :    USE input_constants,                 ONLY: do_gapw_log
      18             :    USE kinds,                           ONLY: dp
      19             :    USE mathconstants,                   ONLY: dfac,&
      20             :                                               twopi
      21             :    USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
      22             :                                               create_grid_atom,&
      23             :                                               grid_atom_type
      24             : #include "./base/base_uses.f90"
      25             : 
      26             :    IMPLICIT NONE
      27             : 
      28             :    PRIVATE
      29             : 
      30             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_set_basis'
      31             : 
      32             :    INTEGER, PARAMETER                                 :: nua = 40, nup = 20
      33             :    REAL(KIND=dp), DIMENSION(nua), PARAMETER :: ugbs = (/0.007299_dp, 0.013705_dp, 0.025733_dp, &
      34             :                                         0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp, &
      35             :                                                  3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
      36             :                                                       174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
      37             :                                                 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
      38             :                                          94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
      39             :                                                       2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
      40             :                                                   27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
      41             :                                                         341890134.751331_dp/)
      42             : 
      43             :    PUBLIC :: set_kind_basis_atomic
      44             : 
      45             : ! **************************************************************************************************
      46             : 
      47             : CONTAINS
      48             : 
      49             : ! **************************************************************************************************
      50             : !> \brief ...
      51             : !> \param basis ...
      52             : !> \param orb_basis_set ...
      53             : !> \param has_pp ...
      54             : !> \param agrid ...
      55             : !> \param cp2k_norm ...
      56             : ! **************************************************************************************************
      57        8458 :    SUBROUTINE set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid, cp2k_norm)
      58             :       TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      59             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
      60             :       LOGICAL, INTENT(IN)                                :: has_pp
      61             :       TYPE(grid_atom_type), OPTIONAL                     :: agrid
      62             :       LOGICAL, INTENT(IN), OPTIONAL                      :: cp2k_norm
      63             : 
      64             :       INTEGER                                            :: i, ii, ipgf, j, k, l, m, ngp, nj, nr, &
      65             :                                                             ns, nset, nsgf, nu, quadtype
      66        8458 :       INTEGER, DIMENSION(:), POINTER                     :: lmax, lmin, npgf, nshell
      67        8458 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf, last_sgf, ls
      68             :       LOGICAL                                            :: has_basis, set_norm
      69             :       REAL(KIND=dp)                                      :: al, an, cn, ear, en, rk
      70        8458 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
      71        8458 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
      72             :       TYPE(grid_atom_type), POINTER                      :: grid
      73             : 
      74        8458 :       IF (ASSOCIATED(orb_basis_set)) THEN
      75             :          has_basis = .TRUE.
      76             :       ELSE
      77           2 :          has_basis = .FALSE.
      78             :       END IF
      79             : 
      80        8458 :       IF (PRESENT(cp2k_norm)) THEN
      81           0 :          set_norm = cp2k_norm
      82             :       ELSE
      83             :          set_norm = .FALSE.
      84             :       END IF
      85             : 
      86        8458 :       NULLIFY (grid)
      87        8458 :       IF (PRESENT(agrid)) THEN
      88           2 :          ngp = agrid%nr
      89           2 :          quadtype = agrid%quadrature
      90             :       ELSE
      91        8456 :          ngp = 400
      92        8456 :          quadtype = do_gapw_log
      93             :       END IF
      94        8458 :       CALL allocate_grid_atom(grid)
      95        8458 :       CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
      96        8458 :       grid%nr = ngp
      97        8458 :       basis%grid => grid
      98             : 
      99        8458 :       NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
     100             : 
     101        8458 :       IF (has_basis) THEN
     102             :          ! fill in the basis data structures
     103        8456 :          basis%basis_type = CGTO_BASIS
     104        8456 :          basis%eps_eig = 1.e-12_dp
     105             :          CALL get_gto_basis_set(orb_basis_set, &
     106             :                                 nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, &
     107             :                                 l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
     108        8456 :                                 first_sgf=first_sgf, last_sgf=last_sgf)
     109       59192 :          basis%nprim = 0
     110       59192 :          basis%nbas = 0
     111       25526 :          DO i = 1, nset
     112       40599 :             DO j = lmin(i), MIN(lmax(i), lmat)
     113       40599 :                basis%nprim(j) = basis%nprim(j) + npgf(i)
     114             :             END DO
     115       58828 :             DO j = 1, nshell(i)
     116       33302 :                l = ls(j, i)
     117       50372 :                IF (l <= lmat) THEN
     118       33302 :                   basis%nbas(l) = basis%nbas(l) + 1
     119       33302 :                   k = basis%nbas(l)
     120             :                END IF
     121             :             END DO
     122             :          END DO
     123             : 
     124       59192 :          nj = MAXVAL(basis%nprim)
     125       59192 :          ns = MAXVAL(basis%nbas)
     126       25362 :          ALLOCATE (basis%am(nj, 0:lmat))
     127      300236 :          basis%am = 0._dp
     128       42268 :          ALLOCATE (basis%cm(nj, ns, 0:lmat))
     129      683924 :          basis%cm = 0._dp
     130       59192 :          DO j = 0, lmat
     131       50736 :             nj = 0
     132       50736 :             ns = 0
     133       50736 :             cn = 2.0_dp**(j + 2)/SQRT(dfac(2*j + 1))/twopi**0.25_dp
     134       50736 :             en = (2*j + 3)*0.25_dp
     135      161612 :             DO i = 1, nset
     136      153156 :                IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
     137       96869 :                   DO ipgf = 1, npgf(i)
     138       96869 :                      basis%am(nj + ipgf, j) = zet(ipgf, i)
     139             :                   END DO
     140       80359 :                   DO ii = 1, nshell(i)
     141       80359 :                      IF (ls(ii, i) == j) THEN
     142       33302 :                         ns = ns + 1
     143       33302 :                         IF (set_norm) THEN
     144           0 :                            DO ipgf = 1, npgf(i)
     145           0 :                               an = cn*zet(ipgf, i)**en
     146           0 :                               basis%cm(nj + ipgf, ns, j) = an*gcc(ipgf, ii, i)
     147             :                            END DO
     148             :                         ELSE
     149      153858 :                            DO ipgf = 1, npgf(i)
     150      153858 :                               basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
     151             :                            END DO
     152             :                         END IF
     153             :                      END IF
     154             :                   END DO
     155       23529 :                   nj = nj + npgf(i)
     156             :                END IF
     157             :             END DO
     158             :          END DO
     159             :          ! Normalization
     160       16912 :          IF (set_norm) THEN
     161           0 :             CALL normalize_basis_cp2k(basis)
     162             :          END IF
     163             :       ELSE
     164             :          ! use default basis
     165           2 :          IF (has_pp) THEN
     166             :             nu = nup
     167             :          ELSE
     168           0 :             nu = nua
     169             :          END IF
     170           2 :          basis%geometrical = .FALSE.
     171           2 :          basis%aval = 0._dp
     172           2 :          basis%cval = 0._dp
     173          14 :          basis%start = 0
     174           2 :          basis%eps_eig = 1.e-12_dp
     175             : 
     176           2 :          basis%basis_type = GTO_BASIS
     177          14 :          basis%nbas = nu
     178          14 :          basis%nprim = nu
     179           4 :          ALLOCATE (basis%am(nu, 0:lmat))
     180          14 :          DO i = 0, lmat
     181         254 :             basis%am(1:nu, i) = ugbs(1:nu)
     182             :          END DO
     183             :       END IF
     184             : 
     185             :       ! initialize basis function on a radial grid
     186        8458 :       nr = basis%grid%nr
     187       59206 :       m = MAXVAL(basis%nbas)
     188       42284 :       ALLOCATE (basis%bf(nr, m, 0:lmat))
     189       25368 :       ALLOCATE (basis%dbf(nr, m, 0:lmat))
     190       25368 :       ALLOCATE (basis%ddbf(nr, m, 0:lmat))
     191             : 
     192    41348014 :       basis%bf = 0._dp
     193    41348014 :       basis%dbf = 0._dp
     194    41348014 :       basis%ddbf = 0._dp
     195       59206 :       DO l = 0, lmat
     196      132786 :          DO i = 1, basis%nprim(l)
     197       73580 :             al = basis%am(i, l)
     198      124328 :             IF (basis%basis_type == GTO_BASIS) THEN
     199      223200 :                DO k = 1, nr
     200      222960 :                   rk = basis%grid%rad(k)
     201      222960 :                   ear = EXP(-al*basis%grid%rad(k)**2)
     202      222960 :                   basis%bf(k, i, l) = rk**l*ear
     203      222960 :                   basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
     204             :                   basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
     205      223200 :                                          2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
     206             :                END DO
     207       73340 :             ELSEIF (basis%basis_type == CGTO_BASIS) THEN
     208    29409340 :                DO k = 1, nr
     209    29336000 :                   rk = basis%grid%rad(k)
     210    29336000 :                   ear = EXP(-al*basis%grid%rad(k)**2)
     211    87630140 :                   DO j = 1, basis%nbas(l)
     212    58220800 :                      basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
     213             :                      basis%dbf(k, j, l) = basis%dbf(k, j, l) &
     214    58220800 :                                           + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
     215             :                      basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
     216             :                                            (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + &
     217    87556800 :                                             4._dp*al*rk**(l + 2))*ear*basis%cm(i, j, l)
     218             :                   END DO
     219             :                END DO
     220             :             ELSE
     221           0 :                CPABORT('Atom basis type?')
     222             :             END IF
     223             :          END DO
     224             :       END DO
     225             : 
     226        8458 :    END SUBROUTINE set_kind_basis_atomic
     227             : 
     228             : ! **************************************************************************************************
     229             : !> \brief ...
     230             : !> \param basis ...
     231             : ! **************************************************************************************************
     232           0 :    SUBROUTINE normalize_basis_cp2k(basis)
     233             :       TYPE(atom_basis_type), INTENT(INOUT)               :: basis
     234             : 
     235             :       INTEGER                                            :: ii, l, n, np
     236             :       REAL(KIND=dp)                                      :: fnorm
     237           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: smat
     238             : 
     239           0 :       DO l = 0, lmat
     240           0 :          n = basis%nbas(l)
     241           0 :          np = basis%nprim(l)
     242           0 :          IF (n > 0) THEN
     243           0 :             ALLOCATE (smat(np, np))
     244           0 :             CALL sg_overlap(smat, l, basis%am(1:np, l), basis%am(1:np, l))
     245           0 :             DO ii = 1, basis%nbas(l)
     246           0 :                fnorm = DOT_PRODUCT(basis%cm(1:np, ii, l), MATMUL(smat, basis%cm(1:np, ii, l)))
     247           0 :                fnorm = 1._dp/SQRT(fnorm)
     248           0 :                basis%cm(1:np, ii, l) = fnorm*basis%cm(1:np, ii, l)
     249             :             END DO
     250           0 :             DEALLOCATE (smat)
     251             :          END IF
     252             :       END DO
     253             : 
     254           0 :    END SUBROUTINE normalize_basis_cp2k
     255             : 
     256             : ! **************************************************************************************************
     257             : 
     258           0 : END MODULE atom_set_basis