LCOV - code coverage report
Current view: top level - src/aobasis - ai_overlap_ppl.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 202 246 82.1 %
Date: 2025-02-18 08:24:35 Functions: 5 6 83.3 %

          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 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 gamma,                           ONLY: fgamma => fgamma_0
      33             :    USE gfun,                            ONLY: gfun_values
      34             :    USE kinds,                           ONLY: dp
      35             :    USE mathconstants,                   ONLY: pi
      36             :    USE mathlib,                         ONLY: binomial
      37             :    USE orbital_pointers,                ONLY: indco,&
      38             :                                               ncoset
      39             : #include "../base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_ppl'
      46             : 
      47             : ! *** Public subroutines ***
      48             : 
      49             :    PUBLIC :: ecploc_integral, ppl_integral, ppl_integral_ri
      50             : 
      51             : CONTAINS
      52             : 
      53             : ! **************************************************************************************************
      54             : !> \brief   Calculation of three-center overlap integrals <a|c|b> over
      55             : !>           Cartesian Gaussian functions for the local part of the Goedecker
      56             : !>           pseudopotential (GTH). c is a primitive Gaussian-type function
      57             : !>           with a set of even angular momentum indices.
      58             : !>
      59             : !>           <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
      60             : !>                               C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
      61             : !>           zetc = alpha**2/2
      62             : !>
      63             : !> \param la_max_set ...
      64             : !> \param la_min_set ...
      65             : !> \param npgfa ...
      66             : !> \param rpgfa ...
      67             : !> \param zeta ...
      68             : !> \param lb_max_set ...
      69             : !> \param lb_min_set ...
      70             : !> \param npgfb ...
      71             : !> \param rpgfb ...
      72             : !> \param zetb ...
      73             : !> \param nexp_ppl ...
      74             : !> \param alpha_ppl ...
      75             : !> \param nct_ppl ...
      76             : !> \param cexp_ppl ...
      77             : !> \param rpgfc ...
      78             : !> \param rab ...
      79             : !> \param dab ...
      80             : !> \param rac ...
      81             : !> \param dac ...
      82             : !> \param rbc ...
      83             : !> \param dbc ...
      84             : !> \param vab ...
      85             : !> \param s ...
      86             : !> \param pab ...
      87             : !> \param force_a ...
      88             : !> \param force_b ...
      89             : !> \param fs ...
      90             : !> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
      91             : !> \param hab2_work ...
      92             : !> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
      93             : !> \param iatom ...
      94             : !> \param jatom ...
      95             : !> \param katom ...
      96             : !> \date    May 2011
      97             : !> \author  Juerg Hutter
      98             : !> \version 1.0
      99             : !> \note    Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
     100             : ! **************************************************************************************************
     101    35009511 :    SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     102    35009511 :                            lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
     103    70019022 :                            rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
     104    70019022 :                            hab2, hab2_work, deltaR, iatom, jatom, katom)
     105             :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     106             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     107             :       INTEGER, INTENT(IN)                                :: lb_max_set, lb_min_set, npgfb
     108             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb, zetb
     109             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     110             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha_ppl
     111             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nct_ppl
     112             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: cexp_ppl
     113             :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     114             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     115             :       REAL(KIND=dp), INTENT(IN)                          :: dab
     116             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     117             :       REAL(KIND=dp), INTENT(IN)                          :: dac
     118             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rbc
     119             :       REAL(KIND=dp), INTENT(IN)                          :: dbc
     120             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     121             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: s
     122             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     123             :          OPTIONAL                                        :: pab
     124             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
     125             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
     126             :          OPTIONAL                                        :: fs, hab2, hab2_work
     127             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     128             :          OPTIONAL                                        :: deltaR
     129             :       INTEGER, INTENT(IN), OPTIONAL                      :: iatom, jatom, katom
     130             : 
     131             :       INTEGER                                            :: iexp, ij, ipgf, jpgf, mmax, nexp
     132             :       REAL(KIND=dp)                                      :: rho, sab, t, zetc
     133    35009511 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: auxint
     134             :       REAL(KIND=dp), DIMENSION(3)                        :: pci
     135             : 
     136    35009511 :       IF (PRESENT(pab)) THEN
     137     9003167 :          CPASSERT(PRESENT(force_a))
     138     9003167 :          CPASSERT(PRESENT(force_b))
     139     9003167 :          CPASSERT(PRESENT(fs))
     140     9003167 :          mmax = la_max_set + lb_max_set + 2
     141     9003167 :          force_a(:) = 0.0_dp
     142     9003167 :          force_b(:) = 0.0_dp
     143    26006344 :       ELSE IF (PRESENT(hab2)) THEN
     144         870 :          mmax = la_max_set + lb_max_set + 2
     145             :       ELSE
     146    26005474 :          mmax = la_max_set + lb_max_set
     147             :       END IF
     148             : 
     149   140038044 :       ALLOCATE (auxint(0:mmax, npgfa*npgfb))
     150  2345194339 :       auxint = 0._dp
     151             : 
     152             :       ! *** Calculate auxiliary integrals ***
     153             : 
     154   155411307 :       DO ipgf = 1, npgfa
     155             :          ! *** Screening ***
     156   120401796 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     157   240038972 :          DO jpgf = 1, npgfb
     158             :             ! *** Screening ***
     159   160214067 :             IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
     160             :                 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
     161    58459095 :             ij = (ipgf - 1)*npgfb + jpgf
     162    58459095 :             rho = zeta(ipgf) + zetb(jpgf)
     163   233836380 :             pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
     164    58459095 :             sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
     165   233836380 :             t = rho*SUM(pci(:)*pci(:))
     166             : 
     167   117182128 :             DO iexp = 1, nexp_ppl
     168    58723033 :                nexp = nct_ppl(iexp)
     169    58723033 :                zetc = alpha_ppl(iexp)
     170   117182128 :                CALL ppl_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
     171             :             END DO
     172             : 
     173   392970023 :             auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
     174             : 
     175             :          END DO
     176             :       END DO
     177             : 
     178             :       CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     179             :                       lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
     180             :                       rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
     181             :                       vab2=hab2, vab2_work=hab2_work, &
     182   192047252 :                       deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
     183             : 
     184    35009511 :       DEALLOCATE (auxint)
     185             : 
     186    35009511 :    END SUBROUTINE ppl_integral
     187             : 
     188             : ! **************************************************************************************************
     189             : !> \brief   Calculation of three-center potential integrals <a|V(r)|b> over
     190             : !>          Cartesian Gaussian functions for the local part of ECP
     191             : !>          pseudopotential.  Multiple terms C1-4 are possible.
     192             : !>
     193             : !>           <a|V(ecploc)|b> = <a| C1/r*exp(-a1*r**2) + C2*exp(-a2*r**2) + C3*r*exp(-a3*r**2) +
     194             : !>                                 C4*r**2*exp(-a4*r**2)|b>
     195             : !>
     196             : !> \param la_max_set ...
     197             : !> \param la_min_set ...
     198             : !> \param npgfa ...
     199             : !> \param rpgfa ...
     200             : !> \param zeta ...
     201             : !> \param lb_max_set ...
     202             : !> \param lb_min_set ...
     203             : !> \param npgfb ...
     204             : !> \param rpgfb ...
     205             : !> \param zetb ...
     206             : !> \param nexp_ppl ...
     207             : !> \param alpha_ppl ...
     208             : !> \param nct_ppl ...
     209             : !> \param cexp_ppl ...
     210             : !> \param rpgfc ...
     211             : !> \param rab ...
     212             : !> \param dab ...
     213             : !> \param rac ...
     214             : !> \param dac ...
     215             : !> \param rbc ...
     216             : !> \param dbc ...
     217             : !> \param vab ...
     218             : !> \param s ...
     219             : !> \param pab ...
     220             : !> \param force_a ...
     221             : !> \param force_b ...
     222             : !> \param fs ...
     223             : !> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
     224             : !> \param hab2_work ...
     225             : !> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
     226             : !> \param iatom ...
     227             : !> \param jatom ...
     228             : !> \param katom ...
     229             : !> \date    2025
     230             : !> \author  Juerg Hutter
     231             : !> \version 1.0
     232             : ! **************************************************************************************************
     233       11380 :    SUBROUTINE ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     234       11380 :                               lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
     235       11380 :                               nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
     236       22760 :                               rab, dab, rac, dac, rbc, dbc, vab, s, pab, &
     237       22760 :                               force_a, force_b, fs, hab2, hab2_work, &
     238       11380 :                               deltaR, iatom, jatom, katom)
     239             :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     240             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     241             :       INTEGER, INTENT(IN)                                :: lb_max_set, lb_min_set, npgfb
     242             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb, zetb
     243             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     244             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha_ppl
     245             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nct_ppl
     246             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: cexp_ppl
     247             :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     248             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     249             :       REAL(KIND=dp), INTENT(IN)                          :: dab
     250             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     251             :       REAL(KIND=dp), INTENT(IN)                          :: dac
     252             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rbc
     253             :       REAL(KIND=dp), INTENT(IN)                          :: dbc
     254             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     255             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: s
     256             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     257             :          OPTIONAL                                        :: pab
     258             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
     259             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
     260             :          OPTIONAL                                        :: fs, hab2, hab2_work
     261             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     262             :          OPTIONAL                                        :: deltaR
     263             :       INTEGER, INTENT(IN), OPTIONAL                      :: iatom, jatom, katom
     264             : 
     265             :       INTEGER                                            :: iexp, ij, ipgf, jpgf, mmax, nexp
     266             :       REAL(KIND=dp)                                      :: rho, sab, t, zetc
     267       11380 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: auxint
     268             :       REAL(KIND=dp), DIMENSION(3)                        :: pci
     269             : 
     270       11380 :       IF (PRESENT(pab)) THEN
     271        2400 :          CPASSERT(PRESENT(force_a))
     272        2400 :          CPASSERT(PRESENT(force_b))
     273        2400 :          CPASSERT(PRESENT(fs))
     274        2400 :          mmax = la_max_set + lb_max_set + 2
     275        2400 :          force_a(:) = 0.0_dp
     276        2400 :          force_b(:) = 0.0_dp
     277        8980 :       ELSE IF (PRESENT(hab2)) THEN
     278           0 :          mmax = la_max_set + lb_max_set + 2
     279             :       ELSE
     280        8980 :          mmax = la_max_set + lb_max_set
     281             :       END IF
     282             : 
     283       45520 :       ALLOCATE (auxint(0:mmax, npgfa*npgfb))
     284      104236 :       auxint = 0._dp
     285             : 
     286             :       ! *** Calculate auxiliary integrals ***
     287             : 
     288       27594 :       DO ipgf = 1, npgfa
     289             :          ! *** Screening ***
     290       16214 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     291       43586 :          DO jpgf = 1, npgfb
     292             :             ! *** Screening ***
     293       19196 :             IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
     294             :                 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
     295       14360 :             ij = (ipgf - 1)*npgfb + jpgf
     296       14360 :             rho = zeta(ipgf) + zetb(jpgf)
     297       57440 :             pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
     298       14360 :             sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
     299       57440 :             t = rho*SUM(pci(:)*pci(:))
     300             : 
     301       34260 :             DO iexp = 1, nexp_ppl
     302       19900 :                nexp = nct_ppl(iexp)
     303       19900 :                zetc = alpha_ppl(iexp)
     304       34260 :                CALL ecploc_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(1, iexp), zetc)
     305             :             END DO
     306             : 
     307       74870 :             auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
     308             : 
     309             :          END DO
     310             :       END DO
     311             : 
     312             :       CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     313             :                       lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
     314             :                       rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
     315             :                       vab2=hab2, vab2_work=hab2_work, &
     316       63480 :                       deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
     317             : 
     318       11380 :       DEALLOCATE (auxint)
     319             : 
     320       11380 :    END SUBROUTINE ecploc_integral
     321             : ! **************************************************************************************************
     322             : !> \brief   Calculation of two-center overlap integrals <a|c> over
     323             : !>          Cartesian Gaussian functions for the local part of the Goedecker
     324             : !>          pseudopotential (GTH). c is a primitive Gaussian-type function
     325             : !>          with a set of even angular momentum indices.
     326             : !>
     327             : !>          <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
     328             : !>                               C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
     329             : !>          zetc = alpha**2/2
     330             : !>
     331             : !> \param la_max_set ...
     332             : !> \param la_min_set ...
     333             : !> \param npgfa ...
     334             : !> \param rpgfa ...
     335             : !> \param zeta ...
     336             : !> \param nexp_ppl ...
     337             : !> \param alpha_ppl ...
     338             : !> \param nct_ppl ...
     339             : !> \param cexp_ppl ...
     340             : !> \param rpgfc ...
     341             : !> \param rac ...
     342             : !> \param dac ...
     343             : !> \param va ...
     344             : !> \param dva ...
     345             : !> \date    December 2017
     346             : !> \author  Juerg Hutter
     347             : !> \version 1.0
     348             : ! **************************************************************************************************
     349         328 :    SUBROUTINE ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     350         328 :                               nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
     351         656 :                               rac, dac, va, dva)
     352             :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     353             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     354             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     355             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha_ppl
     356             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nct_ppl
     357             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: cexp_ppl
     358             :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     359             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     360             :       REAL(KIND=dp), INTENT(IN)                          :: dac
     361             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: va
     362             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
     363             :          OPTIONAL                                        :: dva
     364             : 
     365             :       INTEGER                                            :: i, iexp, ipgf, iw, mmax, na, nexp
     366             :       INTEGER, DIMENSION(3)                              :: ani, anm, anp
     367             :       LOGICAL                                            :: debug
     368             :       REAL(KIND=dp)                                      :: oint, oref, rho, t, zetc
     369         328 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: auxint
     370             :       REAL(KIND=dp), DIMENSION(3)                        :: doint, doref
     371             : 
     372         328 :       debug = .FALSE.
     373             : 
     374         328 :       IF (PRESENT(dva)) THEN
     375         164 :          mmax = la_max_set + 1
     376             :       ELSE
     377         164 :          mmax = la_max_set
     378             :       END IF
     379             : 
     380        1312 :       ALLOCATE (auxint(0:mmax, npgfa))
     381        1588 :       auxint = 0._dp
     382             : 
     383             :       ! *** Calculate auxiliary integrals ***
     384         656 :       DO ipgf = 1, npgfa
     385         328 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     386         328 :          rho = zeta(ipgf)
     387         328 :          t = rho*dac*dac
     388             : 
     389         984 :          DO iexp = 1, nexp_ppl
     390         328 :             nexp = nct_ppl(iexp)
     391         328 :             zetc = alpha_ppl(iexp)
     392         656 :             CALL ppl_aux(auxint(0:mmax, ipgf), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
     393             :          END DO
     394             : 
     395             :       END DO
     396             : 
     397         328 :       IF (PRESENT(dva)) THEN
     398             :          CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     399         164 :                          auxint, rpgfc, rac, dac, va, dva)
     400             :       ELSE
     401             :          CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     402         164 :                          auxint, rpgfc, rac, dac, va)
     403             :       END IF
     404             : 
     405         328 :       DEALLOCATE (auxint)
     406             : 
     407             :       IF (debug) THEN
     408             :          iw = 6
     409             :          na = 0
     410             :          DO ipgf = 1, npgfa
     411             :             IF (rpgfa(ipgf) + rpgfc < dac) THEN
     412             :                na = na + ncoset(la_max_set)
     413             :                CYCLE
     414             :             END IF
     415             :             rho = zeta(ipgf)
     416             :             DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
     417             :                oref = va(na + i)
     418             :                ani(1:3) = indco(1:3, i)
     419             :                oint = ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     420             :                ! test
     421             :                IF (ABS(oint - oref) > 1.0e-12_dp) THEN
     422             :                   WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') "PPL int error     ", ani, la_max_set, dac, oint, oref
     423             :                END IF
     424             :                IF (PRESENT(dva)) THEN
     425             :                   anp = ani + (/1, 0, 0/)
     426             :                   anm = ani - (/1, 0, 0/)
     427             :                   doint(1) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
     428             :                              - ani(1)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     429             :                   anp = ani + (/0, 1, 0/)
     430             :                   anm = ani - (/0, 1, 0/)
     431             :                   doint(2) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
     432             :                              - ani(2)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     433             :                   anp = ani + (/0, 0, 1/)
     434             :                   anm = ani - (/0, 0, 1/)
     435             :                   doint(3) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
     436             :                              - ani(3)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
     437             :                   doref(1:3) = dva(na + i, 1:3)
     438             :                   IF (ANY(ABS(doint - doref) > 1.0e-6_dp)) THEN
     439             :                      WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') " PPL dint error   ", &
     440             :                         ani, la_max_set, dac, SUM(ABS(doint)), SUM(ABS(doref))
     441             :                   END IF
     442             :                END IF
     443             :             END DO
     444             :             na = na + ncoset(la_max_set)
     445             :          END DO
     446             :       END IF
     447             : 
     448         328 :    END SUBROUTINE ppl_integral_ri
     449             : 
     450             : ! **************************************************************************************************
     451             : !> \brief ...
     452             : !> \param rho ...
     453             : !> \param ani ...
     454             : !> \param rac ...
     455             : !> \param nexp_ppl ...
     456             : !> \param nct_ppl ...
     457             : !> \param alpha_ppl ...
     458             : !> \param cexp_ppl ...
     459             : !> \return ...
     460             : ! **************************************************************************************************
     461           0 :    FUNCTION ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) RESULT(oint)
     462             :       REAL(KIND=dp), INTENT(IN)                          :: rho
     463             :       INTEGER, DIMENSION(3), INTENT(IN)                  :: ani
     464             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     465             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     466             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nct_ppl
     467             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha_ppl
     468             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: cexp_ppl
     469             :       REAL(KIND=dp)                                      :: oint
     470             : 
     471             :       INTEGER                                            :: iexp, nexp, ni
     472             :       REAL(KIND=dp)                                      :: cn, zetc
     473             :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     474             : 
     475           0 :       oint = 0.0_dp
     476           0 :       ra = 0.0_dp
     477           0 :       DO iexp = 1, nexp_ppl
     478           0 :          nexp = nct_ppl(iexp)
     479           0 :          zetc = alpha_ppl(iexp)
     480           0 :          CALL init_os_overlap2(rho, zetc, ra, -rac)
     481           0 :          DO ni = 1, nexp
     482           0 :             cn = cexp_ppl(ni, iexp)
     483           0 :             SELECT CASE (ni)
     484             :             CASE (1)
     485           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 0/))
     486             :             CASE (2)
     487           0 :                oint = oint + cn*os_overlap2(ani, (/2, 0, 0/))
     488           0 :                oint = oint + cn*os_overlap2(ani, (/0, 2, 0/))
     489           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 2/))
     490             :             CASE (3)
     491           0 :                oint = oint + cn*os_overlap2(ani, (/4, 0, 0/))
     492           0 :                oint = oint + cn*os_overlap2(ani, (/0, 4, 0/))
     493           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 4/))
     494           0 :                oint = oint + 2.0_dp*cn*os_overlap2(ani, (/2, 2, 0/))
     495           0 :                oint = oint + 2.0_dp*cn*os_overlap2(ani, (/0, 2, 2/))
     496           0 :                oint = oint + 2.0_dp*cn*os_overlap2(ani, (/2, 0, 2/))
     497             :             CASE (4)
     498           0 :                oint = oint + cn*os_overlap2(ani, (/6, 0, 0/))
     499           0 :                oint = oint + cn*os_overlap2(ani, (/0, 6, 0/))
     500           0 :                oint = oint + cn*os_overlap2(ani, (/0, 0, 6/))
     501           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/4, 2, 0/))
     502           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/4, 0, 2/))
     503           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/2, 4, 0/))
     504           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/0, 4, 2/))
     505           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/2, 0, 4/))
     506           0 :                oint = oint + 3.0_dp*cn*os_overlap2(ani, (/0, 2, 4/))
     507           0 :                oint = oint + 6.0_dp*cn*os_overlap2(ani, (/2, 2, 2/))
     508             :             CASE DEFAULT
     509           0 :                CPABORT("OVERLAP_PPL")
     510             :             END SELECT
     511             :          END DO
     512             :       END DO
     513             : 
     514           0 :    END FUNCTION ppl_ri_test
     515             : 
     516             : ! **************************************************************************************************
     517             : !> \brief ...
     518             : !> \param auxint ...
     519             : !> \param mmax ...
     520             : !> \param t ...
     521             : !> \param rho ...
     522             : !> \param nexp_ppl ...
     523             : !> \param cexp_ppl ...
     524             : !> \param zetc ...
     525             : ! **************************************************************************************************
     526    58723361 :    SUBROUTINE ppl_aux(auxint, mmax, t, rho, nexp_ppl, cexp_ppl, zetc)
     527             :       INTEGER, INTENT(IN)                                :: mmax
     528             :       REAL(KIND=dp), DIMENSION(0:mmax)                   :: auxint
     529             :       REAL(KIND=dp), INTENT(IN)                          :: t, rho
     530             :       INTEGER, INTENT(IN)                                :: nexp_ppl
     531             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: cexp_ppl
     532             :       REAL(KIND=dp), INTENT(IN)                          :: zetc
     533             : 
     534             :       INTEGER                                            :: i, j, ke, kp, pmax
     535             :       REAL(KIND=dp)                                      :: a2, a3, a4, cc, f, q, q2, q4, q6, rho2, &
     536             :                                                             rho3, t2, t3
     537             :       REAL(KIND=dp), DIMENSION(0:6)                      :: polder
     538   117446722 :       REAL(KIND=dp), DIMENSION(0:mmax)                   :: expder
     539             : 
     540    58723361 :       CPASSERT(nexp_ppl > 0)
     541    58723361 :       q = rho + zetc
     542    58723361 :       polder = 0._dp
     543    58723361 :       pmax = 0
     544    58723361 :       IF (nexp_ppl > 0) THEN
     545    58723361 :          polder(0) = polder(0) + cexp_ppl(1)
     546    58723361 :          pmax = 0
     547             :       END IF
     548    58723361 :       IF (nexp_ppl > 1) THEN
     549    44453308 :          q2 = q*q
     550    44453308 :          a2 = 0.5_dp/q2*cexp_ppl(2)
     551    44453308 :          polder(0) = polder(0) + a2*(2._dp*rho*t + 3._dp*q)
     552    44453308 :          polder(1) = polder(1) - a2*2._dp*rho
     553    44453308 :          pmax = 1
     554             :       END IF
     555    58723361 :       IF (nexp_ppl > 2) THEN
     556     1211770 :          q4 = q2*q2
     557     1211770 :          rho2 = rho*rho
     558     1211770 :          t2 = t*t
     559     1211770 :          a3 = 0.25_dp/q4*cexp_ppl(3)
     560     1211770 :          polder(0) = polder(0) + a3*(4._dp*rho2*t2 + 20._dp*rho*t*q + 15._dp*q2)
     561     1211770 :          polder(1) = polder(1) - a3*(8._dp*rho2*t + 20._dp*rho*q)
     562     1211770 :          polder(2) = polder(2) + a3*8._dp*rho2
     563     1211770 :          pmax = 2
     564             :       END IF
     565    58723361 :       IF (nexp_ppl > 3) THEN
     566     1211770 :          q6 = q4*q2
     567     1211770 :          rho3 = rho2*rho
     568     1211770 :          t3 = t2*t
     569     1211770 :          a4 = 0.125_dp/q6*cexp_ppl(4)
     570     1211770 :          polder(0) = polder(0) + a4*(8._dp*rho3*t3 + 84._dp*rho2*t2*q + 210._dp*rho*t*q2 + 105._dp*q*q2)
     571     1211770 :          polder(1) = polder(1) - a4*(24._dp*rho3*t2 + 168._dp*rho2*t*q + 210._dp*rho*q2)
     572     1211770 :          polder(2) = polder(2) + a4*(48._dp*rho3*t + 168._dp*rho2*q)
     573     1211770 :          polder(3) = polder(3) - a4*48_dp*rho3
     574     1211770 :          pmax = 3
     575             :       END IF
     576    58723361 :       IF (nexp_ppl > 4) THEN
     577           0 :          CPABORT("nexp_ppl > 4")
     578             :       END IF
     579             : 
     580    58723361 :       f = zetc/q
     581    58723361 :       cc = (pi/q)**1.5_dp*EXP(-t*f)
     582             : 
     583    58723361 :       IF (mmax >= 0) expder(0) = cc
     584   213844661 :       DO i = 1, mmax
     585   213844661 :          expder(i) = f*expder(i - 1)
     586             :       END DO
     587             : 
     588   272568022 :       DO i = 0, mmax
     589   603708685 :          DO j = 0, MIN(i, pmax)
     590   331140663 :             kp = j
     591   331140663 :             ke = i - j
     592   544985324 :             auxint(i) = auxint(i) + expder(ke)*polder(kp)*binomial(i, j)
     593             :          END DO
     594             :       END DO
     595             : 
     596    58723361 :    END SUBROUTINE ppl_aux
     597             : ! **************************************************************************************************
     598             : !> \brief ...
     599             : !> \param auxint ...
     600             : !> \param mmax ...
     601             : !> \param t ...
     602             : !> \param rho ...
     603             : !> \param nexp ...
     604             : !> \param cexp ...
     605             : !> \param zetc ...
     606             : ! **************************************************************************************************
     607       19900 :    SUBROUTINE ecploc_aux(auxint, mmax, t, rho, nexp, cexp, zetc)
     608             :       INTEGER, INTENT(IN)                                :: mmax
     609             :       REAL(KIND=dp), DIMENSION(0:mmax)                   :: auxint
     610             :       REAL(KIND=dp), INTENT(IN)                          :: t, rho
     611             :       INTEGER, INTENT(IN)                                :: nexp
     612             :       REAL(KIND=dp), INTENT(IN)                          :: cexp, zetc
     613             : 
     614             :       INTEGER                                            :: i, j, ke, kf
     615             :       REAL(KIND=dp)                                      :: c0, c1, cc, cval, fa, fr, q, ts
     616       19900 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: expder, fdiff, funder, gfund
     617             : 
     618       19900 :       q = rho + zetc
     619       19900 :       fa = zetc/q
     620       19900 :       fr = rho/q
     621             :       !
     622       99500 :       ALLOCATE (expder(0:mmax), funder(0:mmax + 1))
     623             :       !
     624       20440 :       SELECT CASE (nexp)
     625             :       CASE (0)
     626         540 :          cval = 2.0_dp*cexp/SQRT(q)*pi**1.5_dp*EXP(-t*fa)
     627         540 :          expder(0) = cval
     628        1296 :          DO i = 1, mmax
     629        1296 :             expder(i) = fa*expder(i - 1)
     630             :          END DO
     631         540 :          ts = fr*t
     632        1080 :          ALLOCATE (gfund(0:mmax))
     633         540 :          CALL gfun_values(mmax, ts, gfund)
     634             : 
     635         540 :          funder(0) = gfund(0)
     636        1296 :          DO i = 1, mmax
     637         756 :             funder(i) = 0.0_dp
     638        3267 :             DO j = 0, i
     639        2727 :                funder(i) = funder(i) + (-1)**j*binomial(i, j)*gfund(j)
     640             :             END DO
     641             :          END DO
     642             : 
     643         540 :          DEALLOCATE (gfund)
     644        1296 :          DO i = 1, mmax
     645        1296 :             funder(i) = fr**i*funder(i)
     646             :          END DO
     647        1836 :          DO i = 0, mmax
     648        4347 :             DO j = 0, i
     649        2511 :                kf = j
     650        2511 :                ke = i - j
     651        3807 :                auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
     652             :             END DO
     653             :          END DO
     654             :       CASE (1)
     655        1690 :          cval = cexp*2._dp*pi/q*EXP(-t*fa)
     656        1690 :          expder(0) = cval
     657        4286 :          DO i = 1, mmax
     658        4286 :             expder(i) = fa*expder(i - 1)
     659             :          END DO
     660        1690 :          ts = fr*t
     661        1690 :          CALL fgamma(mmax, ts, funder)
     662        4286 :          DO i = 1, mmax
     663        4286 :             funder(i) = fr**i*funder(i)
     664             :          END DO
     665        5976 :          DO i = 0, mmax
     666       14734 :             DO j = 0, i
     667        8758 :                kf = j
     668        8758 :                ke = i - j
     669       13044 :                auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
     670             :             END DO
     671             :          END DO
     672             :       CASE (2)
     673       17060 :          cval = cexp*(pi/q)**1.5_dp*EXP(-t*fa)
     674       17060 :          expder(0) = cval
     675       49608 :          DO i = 1, mmax
     676       49608 :             expder(i) = fa*expder(i - 1)
     677             :          END DO
     678       66668 :          auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
     679             :       CASE (3)
     680         610 :          cval = 2.*pi*cexp/q**2*EXP(-t*fa)
     681         610 :          expder(0) = cval
     682        1694 :          DO i = 1, mmax
     683        1694 :             expder(i) = fa*expder(i - 1)
     684             :          END DO
     685         610 :          ts = fr*t
     686         610 :          CALL fgamma(mmax + 1, ts, funder)
     687        1220 :          ALLOCATE (fdiff(0:mmax))
     688         610 :          fdiff(0) = (1.0_dp + ts)*funder(0) - ts*funder(1)
     689        1694 :          DO i = 1, mmax
     690             :             fdiff(i) = fr**i*(-i*funder(i - 1) + (1.0_dp + ts)*funder(i) &
     691        1694 :                               + i*funder(i) - ts*funder(i + 1))
     692             :          END DO
     693        2304 :          DO i = 0, mmax
     694        6040 :             DO j = 0, i
     695        3736 :                kf = j
     696        3736 :                ke = i - j
     697        5430 :                auxint(i) = auxint(i) + expder(ke)*fdiff(kf)*binomial(i, j)
     698             :             END DO
     699             :          END DO
     700         610 :          DEALLOCATE (fdiff)
     701             :       CASE (4)
     702           0 :          cval = cexp/(4._dp*q**2)*(pi/q)**1.5_dp*EXP(-t*fa)
     703           0 :          expder(0) = cval
     704           0 :          DO i = 1, mmax
     705           0 :             expder(i) = fa*expder(i - 1)
     706             :          END DO
     707           0 :          c0 = 4._dp*rho/fa
     708           0 :          c1 = 6._dp*q + 4._dp*rho*t
     709           0 :          DO i = 0, mmax
     710           0 :             cc = -i*c0 + c1
     711           0 :             expder(i) = cc*expder(i)
     712             :          END DO
     713           0 :          auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
     714             :       CASE DEFAULT
     715       19900 :          CPABORT("nexp out of range [1..4]")
     716             :       END SELECT
     717             :       !
     718       19900 :       DEALLOCATE (expder, funder)
     719             : 
     720       19900 :    END SUBROUTINE ecploc_aux
     721             : ! **************************************************************************************************
     722             : 
     723             : END MODULE ai_overlap_ppl

Generated by: LCOV version 1.15