Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Calculation of three-center overlap integrals over Cartesian
      10             : !>      Gaussian-type functions for the second term V(ppl) of the local
      11             : !>      part of the Goedecker pseudopotential (GTH):
      12             : !>
      13             : !>      <a|V(local)|b> = <a|V(erf) + V(ppl)|b>
      14             : !>                     = <a|V(erf)|b> + <a|V(ppl)|b>
      15             : !>                     = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
      16             : !>                       (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
      17             : !>                        C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
      18             : !> \par Literature
      19             : !>      S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
      20             : !>      S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
      21             : !>      C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
      22             : !> \par History
      23             : !>      - Derivatives added (17.05.2002,MK)
      24             : !>      - Complete refactoring (05.2011,jhu)
      25             : !> \author Matthias Krack (04.10.2000)
      26             : ! **************************************************************************************************
      27             : MODULE ai_overlap_ppl
      28             :    USE ai_oneelectron,                  ONLY: os_2center,&
      29             :                                               os_3center
      30             :    USE ai_overlap_debug,                ONLY: init_os_overlap2,&
      31             :                                               os_overlap2
      32             :    USE kinds,                           ONLY: dp
      33             :    USE mathconstants,                   ONLY: fac,&
      34             :                                               pi
      35             :    USE orbital_pointers,                ONLY: indco,&
      36             :                                               ncoset
      37             : #include "../base/base_uses.f90"
      38             : 
      39             :    IMPLICIT NONE
      40             : 
      41             :    PRIVATE
      42             : 
      43             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_ppl'
      44             : 
      45             : ! *** Public subroutines ***
      46             : 
      47             :    PUBLIC :: ppl_integral, ppl_integral_ri
      48             : 
      49             : CONTAINS
      50             : 
      51             : ! **************************************************************************************************
      52             : !> \brief   Calculation of three-center overlap integrals <a|c|b> over
      53             : !>           Cartesian Gaussian functions for the local part of the Goedecker
      54             : !>           pseudopotential (GTH). c is a primitive Gaussian-type function
      55             : !>           with a set of even angular momentum indices.
      56             : !>
      57             : !>           <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
      58             : !>                               C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
      59             : !>           zetc = alpha**2/2
      60             : !>
      61             : !> \param la_max_set ...
      62             : !> \param la_min_set ...
      63             : !> \param npgfa ...
      64             : !> \param rpgfa ...
      65             : !> \param zeta ...
      66             : !> \param lb_max_set ...
      67             : !> \param lb_min_set ...
      68             : !> \param npgfb ...
      69             : !> \param rpgfb ...
      70             : !> \param zetb ...
      71             : !> \param nexp_ppl ...
      72             : !> \param alpha_ppl ...
      73             : !> \param nct_ppl ...
      74             : !> \param cexp_ppl ...
      75             : !> \param rpgfc ...
      76             : !> \param rab ...
      77             : !> \param dab ...
      78             : !> \param rac ...
      79             : !> \param dac ...
      80             : !> \param rbc ...
      81             : !> \param dbc ...
      82             : !> \param vab ...
      83             : !> \param s ...
      84             : !> \param pab ...
      85             : !> \param force_a ...
      86             : !> \param force_b ...
      87             : !> \param fs ...
      88             : !> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
      89             : !> \param hab2_work ...
      90             : !> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
      91             : !> \param iatom ...
      92             : !> \param jatom ...
      93             : !> \param katom ...
      94             : !> \date    May 2011
      95             : !> \author  Juerg Hutter
      96             : !> \version 1.0
      97             : !> \note    Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
      98             : ! **************************************************************************************************
      99    34492497 :    SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     100    34492497 :                            lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
     101    68984994 :                            rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
     102    68984994 :                            hab2, hab2_work, deltaR, iatom, jatom, katom)
     103             :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     104             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     105             :       INTEGER, INTENT(IN)                                :: lb_max_set, lb_min_set, npgfb
     106             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb, zetb
     107             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     108             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha_ppl
     109             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nct_ppl
     110             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: cexp_ppl
     111             :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     112             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     113             :       REAL(KIND=dp), INTENT(IN)                          :: dab
     114             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     115             :       REAL(KIND=dp), INTENT(IN)                          :: dac
     116             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rbc
     117             :       REAL(KIND=dp), INTENT(IN)                          :: dbc
     118             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     119             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: s
     120             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     121             :          OPTIONAL                                        :: pab
     122             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
     123             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
     124             :          OPTIONAL                                        :: fs, hab2, hab2_work
     125             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     126             :          OPTIONAL                                        :: deltaR
     127             :       INTEGER, INTENT(IN), OPTIONAL                      :: iatom, jatom, katom
     128             : 
     129             :       INTEGER                                            :: iexp, ij, ipgf, jpgf, mmax, nexp
     130             :       REAL(KIND=dp)                                      :: rho, sab, t, zetc
     131    34492497 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: auxint
     132             :       REAL(KIND=dp), DIMENSION(3)                        :: pci
     133             : 
     134    34492497 :       IF (PRESENT(pab)) THEN
     135     8774386 :          CPASSERT(PRESENT(force_a))
     136     8774386 :          CPASSERT(PRESENT(force_b))
     137     8774386 :          CPASSERT(PRESENT(fs))
     138     8774386 :          mmax = la_max_set + lb_max_set + 2
     139     8774386 :          force_a(:) = 0.0_dp
     140     8774386 :          force_b(:) = 0.0_dp
     141    25718111 :       ELSE IF (PRESENT(hab2)) THEN
     142         870 :          mmax = la_max_set + lb_max_set + 2
     143             :       ELSE
     144    25717241 :          mmax = la_max_set + lb_max_set
     145             :       END IF
     146             : 
     147   137969988 :       ALLOCATE (auxint(0:mmax, npgfa*npgfb))
     148  2290447636 :       auxint = 0._dp
     149             : 
     150             :       ! *** Calculate auxiliary integrals ***
     151             : 
     152   152671131 :       DO ipgf = 1, npgfa
     153             :          ! *** Screening ***
     154   118178634 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     155   235947794 :          DO jpgf = 1, npgfb
     156             :             ! *** Screening ***
     157   157298509 :             IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
     158             :                 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
     159    57615769 :             ij = (ipgf - 1)*npgfb + jpgf
     160    57615769 :             rho = zeta(ipgf) + zetb(jpgf)
     161   230463076 :             pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
     162    57615769 :             sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
     163   230463076 :             t = rho*SUM(pci(:)*pci(:))
     164             : 
     165   115495476 :             DO iexp = 1, nexp_ppl
     166    57879707 :                nexp = nct_ppl(iexp)
     167    57879707 :                zetc = alpha_ppl(iexp)
     168   115495476 :                CALL ppl_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
     169             :             END DO
     170             : 
     171   385989933 :             auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
     172             : 
     173             :          END DO
     174             :       END DO
     175             : 
     176             :       CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     177             :                       lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
     178             :                       rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
     179             :                       vab2=hab2, vab2_work=hab2_work, &
     180   189402730 :                       deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
     181             : 
     182    34492497 :       DEALLOCATE (auxint)
     183             : 
     184    34492497 :    END SUBROUTINE ppl_integral
     185             : ! **************************************************************************************************
     186             : !> \brief   Calculation of two-center overlap integrals <a|c> over
     187             : !>          Cartesian Gaussian functions for the local part of the Goedecker
     188             : !>          pseudopotential (GTH). c is a primitive Gaussian-type function
     189             : !>          with a set of even angular momentum indices.
     190             : !>
     191             : !>          <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
     192             : !>                               C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
     193             : !>          zetc = alpha**2/2
     194             : !>
     195             : !> \param la_max_set ...
     196             : !> \param la_min_set ...
     197             : !> \param npgfa ...
     198             : !> \param rpgfa ...
     199             : !> \param zeta ...
     200             : !> \param nexp_ppl ...
     201             : !> \param alpha_ppl ...
     202             : !> \param nct_ppl ...
     203             : !> \param cexp_ppl ...
     204             : !> \param rpgfc ...
     205             : !> \param rac ...
     206             : !> \param dac ...
     207             : !> \param va ...
     208             : !> \param dva ...
     209             : !> \date    December 2017
     210             : !> \author  Juerg Hutter
     211             : !> \version 1.0
     212             : ! **************************************************************************************************
     213         328 :    SUBROUTINE ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     214         328 :                               nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
     215         656 :                               rac, dac, va, dva)
     216             :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     217             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     218             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     219             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha_ppl
     220             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nct_ppl
     221             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: cexp_ppl
     222             :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     223             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     224             :       REAL(KIND=dp), INTENT(IN)                          :: dac
     225             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: va
     226             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
     227             :          OPTIONAL                                        :: dva
     228             : 
     229             :       INTEGER                                            :: i, iexp, ipgf, iw, mmax, na, nexp
     230             :       INTEGER, DIMENSION(3)                              :: ani, anm, anp
     231             :       LOGICAL                                            :: debug
     232             :       REAL(KIND=dp)                                      :: oint, oref, rho, t, zetc
     233         328 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: auxint
     234             :       REAL(KIND=dp), DIMENSION(3)                        :: doint, doref
     235             : 
     236         328 :       debug = .FALSE.
     237             : 
     238         328 :       IF (PRESENT(dva)) THEN
     239         164 :          mmax = la_max_set + 1
     240             :       ELSE
     241         164 :          mmax = la_max_set
     242             :       END IF
     243             : 
     244        1312 :       ALLOCATE (auxint(0:mmax, npgfa))
     245        1588 :       auxint = 0._dp
     246             : 
     247             :       ! *** Calculate auxiliary integrals ***
     248         656 :       DO ipgf = 1, npgfa
     249         328 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     250         328 :          rho = zeta(ipgf)
     251         328 :          t = rho*dac*dac
     252             : 
     253         984 :          DO iexp = 1, nexp_ppl
     254         328 :             nexp = nct_ppl(iexp)
     255         328 :             zetc = alpha_ppl(iexp)
     256         656 :             CALL ppl_aux(auxint(0:mmax, ipgf), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
     257             :          END DO
     258             : 
     259             :       END DO
     260             : 
     261         328 :       IF (PRESENT(dva)) THEN
     262             :          CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     263         164 :                          auxint, rpgfc, rac, dac, va, dva)
     264             :       ELSE
     265             :          CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     266         164 :                          auxint, rpgfc, rac, dac, va)
     267             :       END IF
     268             : 
     269         328 :       DEALLOCATE (auxint)
     270             : 
     271             :       IF (debug) THEN
     272             :          iw = 6
     273             :          na = 0
     274             :          DO ipgf = 1, npgfa
     275             :             IF (rpgfa(ipgf) + rpgfc < dac) THEN
     276             :                na = na + ncoset(la_max_set)
     277             :                CYCLE
     278             :             END IF
     279             :             rho = zeta(ipgf)
     280             :             DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
     281             :                oref = va(na + i)
     282             :                ani(1:3) = indco(1:3, i)
     283             :                oint = ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     284             :                ! test
     285             :                IF (ABS(oint - oref) > 1.0e-12_dp) THEN
     286             :                   WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') "PPL int error     ", ani, la_max_set, dac, oint, oref
     287             :                END IF
     288             :                IF (PRESENT(dva)) THEN
     289             :                   anp = ani + (/1, 0, 0/)
     290             :                   anm = ani - (/1, 0, 0/)
     291             :                   doint(1) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
     292             :                              - ani(1)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     293             :                   anp = ani + (/0, 1, 0/)
     294             :                   anm = ani - (/0, 1, 0/)
     295             :                   doint(2) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
     296             :                              - ani(2)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     297             :                   anp = ani + (/0, 0, 1/)
     298             :                   anm = ani - (/0, 0, 1/)
     299             :                   doint(3) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
     300             :                              - ani(3)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     301             :                   doref(1:3) = dva(na + i, 1:3)
     302             :                   IF (ANY(ABS(doint - doref) > 1.0e-6_dp)) THEN
     303             :                      WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') " PPL dint error   ", &
     304             :                         ani, la_max_set, dac, SUM(ABS(doint)), SUM(ABS(doref))
     305             :                   END IF
     306             :                END IF
     307             :             END DO
     308             :             na = na + ncoset(la_max_set)
     309             :          END DO
     310             :       END IF
     311             : 
     312         328 :    END SUBROUTINE ppl_integral_ri
     313             : 
     314             : ! **************************************************************************************************
     315             : !> \brief ...
     316             : !> \param rho ...
     317             : !> \param ani ...
     318             : !> \param rac ...
     319             : !> \param nexp_ppl ...
     320             : !> \param nct_ppl ...
     321             : !> \param alpha_ppl ...
     322             : !> \param cexp_ppl ...
     323             : !> \return ...
     324             : ! **************************************************************************************************
     325           0 :    FUNCTION ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) RESULT(oint)
     326             :       REAL(KIND=dp), INTENT(IN)                          :: rho
     327             :       INTEGER, DIMENSION(3), INTENT(IN)                  :: ani
     328             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     329             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     330             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nct_ppl
     331             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha_ppl
     332             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: cexp_ppl
     333             :       REAL(KIND=dp)                                      :: oint
     334             : 
     335             :       INTEGER                                            :: iexp, nexp, ni
     336             :       REAL(KIND=dp)                                      :: cn, zetc
     337             :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     338             : 
     339           0 :       oint = 0.0_dp
     340           0 :       ra = 0.0_dp
     341           0 :       DO iexp = 1, nexp_ppl
     342           0 :          nexp = nct_ppl(iexp)
     343           0 :          zetc = alpha_ppl(iexp)
     344           0 :          CALL init_os_overlap2(rho, zetc, ra, -rac)
     345           0 :          DO ni = 1, nexp
     346           0 :             cn = cexp_ppl(ni, iexp)
     347           0 :             SELECT CASE (ni)
     348             :             CASE (1)
     349           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 0/))
     350             :             CASE (2)
     351           0 :                oint = oint + cn*os_overlap2(ani, (/2, 0, 0/))
     352           0 :                oint = oint + cn*os_overlap2(ani, (/0, 2, 0/))
     353           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 2/))
     354             :             CASE (3)
     355           0 :                oint = oint + cn*os_overlap2(ani, (/4, 0, 0/))
     356           0 :                oint = oint + cn*os_overlap2(ani, (/0, 4, 0/))
     357           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 4/))
     358           0 :                oint = oint + 2.0_dp*cn*os_overlap2(ani, (/2, 2, 0/))
     359           0 :                oint = oint + 2.0_dp*cn*os_overlap2(ani, (/0, 2, 2/))
     360           0 :                oint = oint + 2.0_dp*cn*os_overlap2(ani, (/2, 0, 2/))
     361             :             CASE (4)
     362           0 :                oint = oint + cn*os_overlap2(ani, (/6, 0, 0/))
     363           0 :                oint = oint + cn*os_overlap2(ani, (/0, 6, 0/))
     364           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 6/))
     365           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/4, 2, 0/))
     366           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/4, 0, 2/))
     367           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/2, 4, 0/))
     368           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/0, 4, 2/))
     369           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/2, 0, 4/))
     370           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/0, 2, 4/))
     371           0 :                oint = oint + 6.0_dp*cn*os_overlap2(ani, (/2, 2, 2/))
     372             :             CASE DEFAULT
     373           0 :                CPABORT("OVERLAP_PPL")
     374             :             END SELECT
     375             :          END DO
     376             :       END DO
     377             : 
     378           0 :    END FUNCTION ppl_ri_test
     379             : 
     380             : ! **************************************************************************************************
     381             : !> \brief ...
     382             : !> \param auxint ...
     383             : !> \param mmax ...
     384             : !> \param t ...
     385             : !> \param rho ...
     386             : !> \param nexp_ppl ...
     387             : !> \param cexp_ppl ...
     388             : !> \param zetc ...
     389             : ! **************************************************************************************************
     390    57880035 :    SUBROUTINE ppl_aux(auxint, mmax, t, rho, nexp_ppl, cexp_ppl, zetc)
     391             :       INTEGER, INTENT(IN)                                :: mmax
     392             :       REAL(KIND=dp), DIMENSION(0:mmax)                   :: auxint
     393             :       REAL(KIND=dp), INTENT(IN)                          :: t, rho
     394             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     395             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: cexp_ppl
     396             :       REAL(KIND=dp), INTENT(IN)                          :: zetc
     397             : 
     398             :       INTEGER                                            :: i, j, ke, kp, pmax
     399             :       REAL(KIND=dp)                                      :: a2, a3, a4, cc, f, q, q2, q4, q6, rho2, &
     400             :                                                             rho3, t2, t3
     401             :       REAL(KIND=dp), DIMENSION(0:6)                      :: polder
     402   115760070 :       REAL(KIND=dp), DIMENSION(0:mmax)                   :: expder
     403             : 
     404    57880035 :       CPASSERT(nexp_ppl > 0)
     405    57880035 :       q = rho + zetc
     406    57880035 :       polder = 0._dp
     407    57880035 :       pmax = 0
     408    57880035 :       IF (nexp_ppl > 0) THEN
     409    57880035 :          polder(0) = polder(0) + cexp_ppl(1)
     410    57880035 :          pmax = 0
     411             :       END IF
     412    57880035 :       IF (nexp_ppl > 1) THEN
     413    43774207 :          q2 = q*q
     414    43774207 :          a2 = 0.5_dp/q2*cexp_ppl(2)
     415    43774207 :          polder(0) = polder(0) + a2*(2._dp*rho*t + 3._dp*q)
     416    43774207 :          polder(1) = polder(1) - a2*2._dp*rho
     417    43774207 :          pmax = 1
     418             :       END IF
     419    57880035 :       IF (nexp_ppl > 2) THEN
     420     1211770 :          q4 = q2*q2
     421     1211770 :          rho2 = rho*rho
     422     1211770 :          t2 = t*t
     423     1211770 :          a3 = 0.25_dp/q4*cexp_ppl(3)
     424     1211770 :          polder(0) = polder(0) + a3*(4._dp*rho2*t2 + 20._dp*rho*t*q + 15._dp*q2)
     425     1211770 :          polder(1) = polder(1) - a3*(8._dp*rho2*t + 20._dp*rho*q)
     426     1211770 :          polder(2) = polder(2) + a3*8._dp*rho2
     427     1211770 :          pmax = 2
     428             :       END IF
     429    57880035 :       IF (nexp_ppl > 3) THEN
     430     1211770 :          q6 = q4*q2
     431     1211770 :          rho3 = rho2*rho
     432     1211770 :          t3 = t2*t
     433     1211770 :          a4 = 0.125_dp/q6*cexp_ppl(4)
     434     1211770 :          polder(0) = polder(0) + a4*(8._dp*rho3*t3 + 84._dp*rho2*t2*q + 210._dp*rho*t*q2 + 105._dp*q*q2)
     435     1211770 :          polder(1) = polder(1) - a4*(24._dp*rho3*t2 + 168._dp*rho2*t*q + 210._dp*rho*q2)
     436     1211770 :          polder(2) = polder(2) + a4*(48._dp*rho3*t + 168._dp*rho2*q)
     437     1211770 :          polder(3) = polder(3) - a4*48_dp*rho3
     438     1211770 :          pmax = 3
     439             :       END IF
     440    57880035 :       IF (nexp_ppl > 4) THEN
     441           0 :          CPABORT("nexp_ppl > 4")
     442             :       END IF
     443             : 
     444    57880035 :       f = zetc/q
     445    57880035 :       cc = (pi/q)**1.5_dp*EXP(-t*f)
     446             : 
     447    57880035 :       IF (mmax >= 0) expder(0) = cc
     448   209940399 :       DO i = 1, mmax
     449   209940399 :          expder(i) = f*expder(i - 1)
     450             :       END DO
     451             : 
     452   267820434 :       DO i = 0, mmax
     453   592819839 :          DO j = 0, MIN(i, pmax)
     454   324999405 :             kp = j
     455   324999405 :             ke = i - j
     456   534939804 :             auxint(i) = auxint(i) + expder(ke)*polder(kp)*choose(i, j)
     457             :          END DO
     458             :       END DO
     459             : 
     460    57880035 :    END SUBROUTINE ppl_aux
     461             : ! **************************************************************************************************
     462             : !> \brief ...
     463             : !> \param n ...
     464             : !> \param k ...
     465             : !> \return ...
     466             : ! **************************************************************************************************
     467   324999405 :    FUNCTION choose(n, k)
     468             : 
     469             :       INTEGER, INTENT(IN)                                :: n, k
     470             :       REAL(KIND=dp)                                      :: choose
     471             : 
     472   324999405 :       IF (n >= k) THEN
     473   324999405 :          choose = REAL(NINT(fac(n)/(fac(k)*fac(n - k))), KIND=dp)
     474             :       ELSE
     475             :          choose = 0.0_dp
     476             :       END IF
     477             : 
     478   324999405 :    END FUNCTION choose
     479             : ! **************************************************************************************************
     480             : 
     481             : END MODULE ai_overlap_ppl