LCOV - code coverage report
Current view: top level - src - pair_potential_util.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.3 % 110 107
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

            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              : !> \par History
      10              : !>      September 2005 - Introduced the Born-Mayer-Huggins-Fumi-Tosi  Potential (BMHTF)
      11              : !>      2006 - Major rewriting of the routines.. Linear scaling setup of splines
      12              : !> \author CJM
      13              : ! **************************************************************************************************
      14              : MODULE pair_potential_util
      15              : 
      16              :    USE fparser,                         ONLY: EvalErrType,&
      17              :                                               evalf
      18              :    USE kinds,                           ONLY: dp
      19              :    USE mathconstants,                   ONLY: ifac
      20              :    USE pair_potential_types,            ONLY: &
      21              :         b4_type, bm_type, ea_type, ft_type, ftd_type, gp_type, gw_type, ip_type, lj_charmm_type, &
      22              :         lj_type, not_initialized, pair_potential_single_type, tab_type, wl_type
      23              :    USE physcon,                         ONLY: bohr,&
      24              :                                               evolt
      25              : #include "./base/base_uses.f90"
      26              : 
      27              :    IMPLICIT NONE
      28              : 
      29              :    PRIVATE
      30              : 
      31              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pair_potential_util'
      32              :    REAL(KIND=dp), PARAMETER, PRIVATE    :: MIN_HICUT_VALUE = 1.0E-15_dp, &
      33              :                                            DEFAULT_HICUT_VALUE = 1.0E3_dp
      34              : 
      35              :    PUBLIC :: ener_pot, ener_zbl, zbl_matching_polinomial
      36              : 
      37              : CONTAINS
      38              : 
      39              : ! **************************************************************************************************
      40              : !> \brief Evaluates the nonbond potential energy for the implemented FF kinds
      41              : !> \param pot ...
      42              : !> \param r ...
      43              : !> \param energy_cutoff ...
      44              : !> \return ...
      45              : ! **************************************************************************************************
      46    779385334 :    FUNCTION ener_pot(pot, r, energy_cutoff) RESULT(value)
      47              :       TYPE(pair_potential_single_type), POINTER          :: pot
      48              :       REAL(KIND=dp), INTENT(IN)                          :: r, energy_cutoff
      49              :       REAL(KIND=dp)                                      :: value
      50              : 
      51              :       INTEGER                                            :: i, index, index1, index2, j, n
      52              :       REAL(KIND=dp)                                      :: bd6, bd8, dampsum, f6, f8, lvalue, pp, &
      53              :                                                             qq, scale, xf
      54              : 
      55    779385334 :       value = 0.0_dp
      56   1558950036 :       DO j = 1, SIZE(pot%type)
      57              :          ! A lower boundary for the potential definition was defined
      58    779564702 :          IF ((pot%set(j)%rmin /= not_initialized) .AND. (r < pot%set(j)%rmin)) CYCLE
      59              :          ! An upper boundary for the potential definition was defined
      60    779564044 :          IF ((pot%set(j)%rmax /= not_initialized) .AND. (r >= pot%set(j)%rmax)) CYCLE
      61              :          ! If within limits let's compute the potential...
      62    779561944 :          IF (pot%type(j) == lj_charmm_type) THEN
      63              :             lvalue = &
      64              :                4.0_dp*pot%set(j)%lj%epsilon*(pot%set(j)%lj%sigma12*r**(-12) - pot%set(j)%lj% &
      65    273070064 :                                              sigma6*r**(-6))
      66              :          ELSE IF (pot%type(j) == lj_type) THEN
      67              :             lvalue = pot%set(j)%lj%epsilon* &
      68       604304 :                      (pot%set(j)%lj%sigma12*r**(-12) - pot%set(j)%lj%sigma6*r**(-6))
      69              :          ELSE IF (pot%type(j) == ip_type) THEN
      70      6057566 :             lvalue = 0._dp
      71      6057566 :             IF (r > pot%set(j)%ipbv%rcore) THEN
      72     66181110 :                DO i = 2, 15
      73     66181110 :                   lvalue = lvalue + pot%set(j)%ipbv%a(i)/(r**(i - 1)*REAL(i - 1, dp))
      74              :                END DO
      75              :             ELSE
      76              :                ! use a linear potential
      77      1645492 :                lvalue = pot%set(j)%ipbv%m*r + pot%set(j)%ipbv%b
      78              :             END IF
      79              :             lvalue = lvalue
      80              :          ELSE IF (pot%type(j) == wl_type) THEN
      81    132970017 :             lvalue = pot%set(j)%willis%a*EXP(-pot%set(j)%willis%b*r) - pot%set(j)%willis%c/r**6
      82              :          ELSE IF (pot%type(j) == gw_type) THEN
      83              :             scale = EXP(pot%set(j)%goodwin%m*(-(r/pot%set(j)%goodwin%dc)**pot%set(j)%goodwin%mc + &
      84            0 :                                               (pot%set(j)%goodwin%d/pot%set(j)%goodwin%dc)**pot%set(j)%goodwin%mc))
      85            0 :             lvalue = scale*pot%set(j)%goodwin%vr0*(pot%set(j)%goodwin%d/r)**pot%set(j)%goodwin%m
      86              :          ELSE IF (pot%type(j) == ft_type) THEN
      87        25762 :             lvalue = pot%set(j)%ft%a*EXP(-pot%set(j)%ft%b*r) - pot%set(j)%ft%c/r**6 - pot%set(j)%ft%d/r**8
      88              :          ELSE IF (pot%type(j) == ftd_type) THEN
      89              :             ! Compute 6th order dispersion correction term
      90      8994460 :             bd6 = pot%set(j)%ftd%bd(1)
      91      8994460 :             dampsum = 1.0_dp
      92      8994460 :             xf = 1.0_dp
      93     62961220 :             DO i = 1, 6
      94     53966760 :                xf = xf*bd6*r
      95     62961220 :                dampsum = dampsum + xf*ifac(i)
      96              :             END DO
      97      8994460 :             f6 = 1.0_dp - EXP(-bd6*r)*dampsum
      98              :             ! Compute 8th order dispersion correction term
      99      8994460 :             bd8 = pot%set(j)%ftd%bd(2)
     100      8994460 :             dampsum = 1.0_dp
     101      8994460 :             xf = 1.0_dp
     102     80950140 :             DO i = 1, 8
     103     71955680 :                xf = xf*bd8*r
     104     80950140 :                dampsum = dampsum + xf*ifac(i)
     105              :             END DO
     106      8994460 :             f8 = 1.0_dp - EXP(-bd8*r)*dampsum
     107      8994460 :             lvalue = pot%set(j)%ftd%a*EXP(-pot%set(j)%ftd%b*r) - f6*pot%set(j)%ftd%c/r**6 - f8*pot%set(j)%ftd%d/r**8
     108              :          ELSE IF (pot%type(j) == ea_type) THEN
     109         3854 :             index = INT(r/pot%set(j)%eam%drar) + 1
     110         3854 :             IF (index > pot%set(j)%eam%npoints) THEN
     111              :                index = pot%set(j)%eam%npoints
     112         3702 :             ELSEIF (index < 1) THEN
     113              :                index = 1
     114              :             END IF
     115         3854 :             qq = r - pot%set(j)%eam%rval(index)
     116              :             pp = pot%set(j)%eam%phi(index) + &
     117         3854 :                  qq*pot%set(j)%eam%phip(index)
     118         3854 :             lvalue = pp
     119              :          ELSE IF (pot%type(j) == b4_type) THEN
     120    285055946 :             IF (r <= pot%set(j)%buck4r%r1) THEN
     121     73446226 :                pp = pot%set(j)%buck4r%a*EXP(-pot%set(j)%buck4r%b*r)
     122    211609720 :             ELSEIF (r > pot%set(j)%buck4r%r1 .AND. r <= pot%set(j)%buck4r%r2) THEN
     123     17904118 :                pp = 0.0_dp
     124    125328826 :                DO n = 0, pot%set(j)%buck4r%npoly1
     125    125328826 :                   pp = pp + pot%set(j)%buck4r%poly1(n)*r**n
     126              :                END DO
     127    193705602 :             ELSEIF (r > pot%set(j)%buck4r%r2 .AND. r <= pot%set(j)%buck4r%r3) THEN
     128     19525932 :                pp = 0.0_dp
     129     97629660 :                DO n = 0, pot%set(j)%buck4r%npoly2
     130     97629660 :                   pp = pp + pot%set(j)%buck4r%poly2(n)*r**n
     131              :                END DO
     132    174179670 :             ELSEIF (r > pot%set(j)%buck4r%r3) THEN
     133    174179670 :                pp = -pot%set(j)%buck4r%c/r**6
     134              :             END IF
     135              :             lvalue = pp
     136              :          ELSE IF (pot%type(j) == tab_type) THEN
     137      5647200 :             index1 = FLOOR((r - pot%set(j)%tab%r(1))/pot%set(j)%tab%dr) + 1
     138      5647200 :             index2 = index1 + 1
     139      5647200 :             IF (index2 > pot%set(j)%tab%npoints) THEN
     140           48 :                index2 = pot%set(j)%tab%npoints
     141           48 :                index1 = index2 - 1
     142      5647152 :             ELSEIF (index1 < 1) THEN
     143       882796 :                index1 = 1
     144       882796 :                index2 = 2
     145              :             END IF
     146              :             pp = pot%set(j)%tab%e(index1) + (r - pot%set(j)%tab%r(index1))* &
     147              :                  (pot%set(j)%tab%e(index2) - pot%set(j)%tab%e(index1))/ &
     148      5647200 :                  (pot%set(j)%tab%r(index2) - pot%set(j)%tab%r(index1))
     149      5647200 :             lvalue = pp
     150              :          ELSE IF (pot%type(j) == bm_type) THEN
     151              :             lvalue = pot%set(j)%buckmo%f0*(pot%set(j)%buckmo%b1 + pot%set(j)%buckmo%b2)* &
     152              :                      EXP((pot%set(j)%buckmo%a1 + pot%set(j)%buckmo%a2 - r)/(pot%set(j)%buckmo%b1 + pot%set(j)%buckmo%b2)) &
     153              :                      - pot%set(j)%buckmo%c/r**6 &
     154              :                      + pot%set(j)%buckmo%d*(EXP(-2._dp*pot%set(j)%buckmo%beta*(r - pot%set(j)%buckmo%r0)) - &
     155       696100 :                                             2.0_dp*EXP(-pot%set(j)%buckmo%beta*(r - pot%set(j)%buckmo%r0)))
     156              :          ELSE IF (pot%type(j) == gp_type) THEN
     157     50359896 :             pot%set(j)%gp%values(1) = r
     158     50359896 :             lvalue = evalf(pot%set(j)%gp%myid, pot%set(j)%gp%values)
     159     50359896 :             IF (EvalErrType > 0) &
     160            0 :                CPABORT("Error evaluating generic potential energy function")
     161              :          ELSE
     162              :             lvalue = 0.0_dp
     163              :          END IF
     164   1558950036 :          value = value + lvalue
     165              :       END DO
     166    779385334 :       value = value - energy_cutoff
     167    779385334 :    END FUNCTION ener_pot
     168              : 
     169              : ! **************************************************************************************************
     170              : !> \brief Evaluates the ZBL scattering potential, very short range
     171              : !>        Only shell-model for interactions among pairs without repulsive term
     172              : !> \param pot ...
     173              : !> \param r ...
     174              : !> \return ...
     175              : !> \author i soliti ignoti
     176              : ! **************************************************************************************************
     177     40236728 :    FUNCTION ener_zbl(pot, r)
     178              : 
     179              :       TYPE(pair_potential_single_type), POINTER          :: pot
     180              :       REAL(KIND=dp), INTENT(IN)                          :: r
     181              :       REAL(KIND=dp)                                      :: ener_zbl
     182              : 
     183              :       REAL(KIND=dp)                                      :: au, fac, x
     184              : 
     185     40236728 :       ener_zbl = 0.0_dp
     186     40236728 :       IF (r <= pot%zbl_rcut(1)) THEN
     187     13579482 :          au = 0.88534_dp*bohr/(pot%z1**0.23_dp + pot%z2**0.23_dp)
     188     13579482 :          x = r/au
     189     13579482 :          fac = pot%z1*pot%z2/evolt
     190              :          ener_zbl = fac/r*(0.1818_dp*EXP(-3.2_dp*x) + 0.5099_dp*EXP(-0.9423_dp*x) + &
     191     13579482 :                            0.2802_dp*EXP(-0.4029_dp*x) + 0.02817_dp*EXP(-0.2016_dp*x))
     192     26657246 :       ELSEIF (r > pot%zbl_rcut(1) .AND. r <= pot%zbl_rcut(2)) THEN
     193              :          ener_zbl = pot%zbl_poly(0) + pot%zbl_poly(1)*r + pot%zbl_poly(2)*r*r + pot%zbl_poly(3)*r*r*r + &
     194      1872516 :                     pot%zbl_poly(4)*r*r*r*r + pot%zbl_poly(5)*r*r*r*r*r
     195              :       ELSE
     196              :          ener_zbl = 0.0_dp
     197              :       END IF
     198              : 
     199     40236728 :    END FUNCTION ener_zbl
     200              : 
     201              : ! **************************************************************************************************
     202              : !> \brief Determine the polinomial coefficients used to set to zero the zbl potential
     203              : !>        at the cutoff radius, with continuity in function, first and second derivative
     204              : !>        Only shell-model for interactions among pairs without repulsive term
     205              : !> \param pot ...
     206              : !> \param rcov1 ...
     207              : !> \param rcov2 ...
     208              : !> \param z1 ...
     209              : !> \param z2 ...
     210              : !> \author i soliti ignoti
     211              : ! **************************************************************************************************
     212           48 :    SUBROUTINE zbl_matching_polinomial(pot, rcov1, rcov2, z1, z2)
     213              : 
     214              :       TYPE(pair_potential_single_type), POINTER          :: pot
     215              :       REAL(KIND=dp), INTENT(IN)                          :: rcov1, rcov2, z1, z2
     216              : 
     217              :       REAL(KIND=dp)                                      :: au, d1, d2, dd1, dd2, fac, v1, v2, x, &
     218              :                                                             x1, x2
     219              : 
     220           48 :       pot%zbl_rcut(1) = (rcov1 + rcov2)*(1.0_dp - 0.2_dp)*bohr
     221           48 :       pot%zbl_rcut(2) = (rcov1 + rcov2)*bohr
     222           48 :       x1 = pot%zbl_rcut(1)
     223           48 :       x2 = pot%zbl_rcut(2)
     224           48 :       pot%z1 = z1
     225           48 :       pot%z2 = z2
     226              : 
     227           48 :       au = 0.88534_dp*bohr/(z1**0.23_dp + z2**0.23_dp)
     228           48 :       x = x1/au
     229           48 :       fac = z1*z2/evolt
     230              :       v1 = fac/x1*(0.1818_dp*EXP(-3.2_dp*x) + 0.5099_dp*EXP(-0.9423_dp*x) + &
     231           48 :                    0.2802_dp*EXP(-0.4029_dp*x) + 0.02817_dp*EXP(-0.2016_dp*x))
     232              :       d1 = fac/x1/au*(-3.2_dp*0.1818_dp*EXP(-3.2_dp*x) - 0.9423_dp*0.5099_dp*EXP(-0.9423_dp*x) &
     233              :                       - 0.4029_dp*0.2802_dp*EXP(-0.4029_dp*x) - 0.2016_dp*0.02817_dp*EXP(-0.2016_dp*x)) &
     234              :            - fac/x1/x1*(0.1818_dp*EXP(-3.2_dp*x) + 0.5099_dp*EXP(-0.9423_dp*x) + &
     235           48 :                         0.2802_dp*EXP(-0.4029_dp*x) + 0.02817_dp*EXP(-0.2016_dp*x))
     236              : 
     237              :       dd1 = 2.0_dp*fac/x1**3*(0.1818_dp*EXP(-0.32E1_dp*x) &
     238              :                               + 0.5099_dp*EXP(-0.9423_dp*x) + 0.2802_dp*EXP(-0.4029_dp*x) &
     239              :                               + 0.2817E-1_dp*EXP(-0.2016_dp*x)) &
     240              :             - 0.2E1_dp*fac/x1**2/au*(-0.58176_dp*EXP(-0.32E1_dp*x) - 0.48047877_dp*EXP(-0.9423_dp*x) &
     241              :                                      - 0.11289258_dp*EXP(-0.4029_dp*x) - 0.5679072E-2_dp*EXP(-0.2016_dp*x)) &
     242              :             + fac/x1/au**2*(0.1861632E1_dp*EXP(-0.32E1_dp*x) + &
     243              :                             0.4527551450_dp*EXP(-0.9423_dp*x) + 0.4548442048E-1_dp*EXP(-0.4029_dp*x) + &
     244           48 :                             0.1144900915E-2_dp*EXP(-0.2016_dp*x))
     245              : 
     246           48 :       v2 = 0.0_dp
     247           48 :       d2 = 0.0_dp
     248           48 :       dd2 = 0.0_dp
     249              : 
     250           48 :       CALL compute_polinomial_5th(x1, v1, d1, dd1, x2, v2, d2, dd2, pot%zbl_poly)
     251              : 
     252           48 :    END SUBROUTINE zbl_matching_polinomial
     253              : 
     254              : ! **************************************************************************************************
     255              : !> \brief ...
     256              : !> \param r1 ...
     257              : !> \param v1 ...
     258              : !> \param d1 ...
     259              : !> \param dd1 ...
     260              : !> \param r2 ...
     261              : !> \param v2 ...
     262              : !> \param d2 ...
     263              : !> \param dd2 ...
     264              : !> \param poly ...
     265              : ! **************************************************************************************************
     266           48 :    SUBROUTINE compute_polinomial_5th(r1, v1, d1, dd1, r2, v2, d2, dd2, poly)
     267              : 
     268              :       REAL(KIND=dp)                                      :: r1, v1, d1, dd1, r2, v2, d2, dd2, &
     269              :                                                             poly(0:5)
     270              : 
     271              :       REAL(KIND=dp)                                      :: a0, a1, a2, a3, a4, a5
     272              : 
     273              : !  5th order
     274              : 
     275              :       a0 = .5_dp*(2._dp*r1**5*v2 - 2._dp*v1*r2**5 + 10._dp*v1*r2**4*r1 - 20._dp*v1*r1**2*r2**3 - r1**2*dd1*r2**5 - &
     276              :                   r1**4*r2**3*dd1 + 20._dp*r1**3*r2**2*v2 + 2._dp*r1**3*r2**4*dd1 + r1**3*r2**4*dd2 - 8._dp*r1**3*r2**3*d2 - &
     277              :                   2._dp*r1**4*r2**3*dd2 + 10._dp*r1**4*r2**2*d2 - 10._dp*r1**4*r2*v2 + 2._dp*r1*d1*r2**5 - 10._dp*r1**2*d1*r2**4 + &
     278              :                   8._dp*r1**3*d1*r2**3 - 2._dp*r1**5*r2*d2 + r1**5*r2**2*dd2)/ &
     279           48 :            (10.*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
     280              : 
     281              :       a1 = -.5_dp*(-4._dp*r2**3*r1**3*dd2 + 24._dp*r2**2*r1**3*d1 + 4._dp*r2**3*r1**3*dd1 + 3._dp*r2**4*r1**2*dd2 + &
     282              :                    r2**4*r1**2*dd1 - 2._dp*r2**5*r1*dd1 - 10._dp*r2**4*r1*d1 + 10._dp*r1**4*r2*d2 - &
     283              :                    r1**4*r2**2*dd2 - 3._dp*r1**4*r2**2*dd1 + &
     284              :                    2._dp*r1**5*r2*dd2 - 24._dp*r2**3*r1**2*d2 - 16._dp*r2**3*r1**2*d1 + &
     285              :                    16._dp*r2**2*r1**3*d2 - 2._dp*r1**5*d2 + 2._dp*r2**5*d1 - &
     286              :                    60._dp*r1**2*r2**2*v1 + 60._dp*r1**2*r2**2*v2)/ &
     287           48 :            (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
     288              : 
     289              :       a2 = .5_dp*(60._dp*r1**2*r2*v2 - 60._dp*v1*r1*r2**2 - 12._dp*r1**2*r2**2*d2 - 36._dp*r1*d1*r2**3 + 3._dp*r2**4*r1*dd2 - &
     290              :                   24._dp*r2**3*r1*d2 - 4._dp*r2**4*r1*dd1 + 12._dp*r1**2*r2**2*d1 - 8._dp*r1**3*r2**2*dd2 + 24._dp*r1**3*r2*d1 + &
     291              :                   4._dp*r1**4*r2*dd2 + 36._dp*r1**3*r2*d2 - 3._dp*r1**4*r2*dd1 + 8._dp*r2**3*r1**2*dd1 + 60._dp*r2**2*r1*v2 - &
     292              :                   60._dp*r1**2*v1*r2 + r1**5*dd2 - r2**5*dd1)/ &
     293           48 :            (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
     294              : 
     295              :       a3 = -.5_dp*(3._dp*r1**4*dd2 - r1**4*dd1 + 8.*r1**3*d1 - 4.*r1**3*r2*dd1 + &
     296              :                    12._dp*r1**3*d2 + 32._dp*r1**2*r2*d1 - 8._dp*r1**2*r2**2*dd2 - &
     297              :                    20._dp*r1**2*v1 + 8._dp*r1**2*r2**2*dd1 + 28._dp*r1**2*r2*d2 + &
     298              :                    20._dp*r1**2*v2 + 80._dp*r1*r2*v2 - 28._dp*r2**2*r1*d1 - 80._dp*r1*v1*r2 - &
     299              :                    32._dp*r2**2*r1*d2 + 4._dp*r1*r2**3*dd2 - 8._dp*r2**3*d2 - 12._dp*r2**3*d1 + &
     300              :                    r2**4*dd2 - 3._dp*r2**4*dd1 + 20._dp*r2**2*v2 - 20._dp*r2**2*v1)/ &
     301           48 :            (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
     302              : 
     303              :       a4 = .5_dp*(3._dp*r1**3*dd2 - 2._dp*r1**3*dd1 + r1**2*r2*dd1 + 14.*r1**2*d1 - 4._dp*r1**2*r2*dd2 + &
     304              :                   16._dp*r1**2*d2 - 2._dp*r1*r2*d2 - r1*r2**2*dd2 - &
     305              :                   30._dp*r1*v1 + 30.*r1*v2 + 2._dp*r1*r2*d1 + 4._dp*r1*r2**2*dd1 - 16._dp*r2**2*d1 + &
     306              :                   2._dp*r2**3*dd2 - 14._dp*r2**2*d2 + 30._dp*r2*v2 - 30._dp*v1*r2 - &
     307              :                   3._dp*r2**3*dd1)/(10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - &
     308           48 :                                     10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
     309              : 
     310              :       a5 = -.5_dp*(6._dp*r1*d1 + 2._dp*r2*r1*dd1 + 6._dp*r1*d2 - 2.*r2*r1*dd2 - &
     311              :                    r2**2*dd1 - r1**2*dd1 - 12.*v1 + 12._dp*v2 + r1**2*dd2 - &
     312              :                    6._dp*r2*d1 + r2**2*dd2 - 6._dp*r2*d2)/ &
     313           48 :            (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
     314              : 
     315           48 :       poly(0) = a0
     316           48 :       poly(1) = a1
     317           48 :       poly(2) = a2
     318           48 :       poly(3) = a3
     319           48 :       poly(4) = a4
     320           48 :       poly(5) = a5
     321              : 
     322           48 :    END SUBROUTINE compute_polinomial_5th
     323              : 
     324              : END MODULE pair_potential_util
     325              : 
        

Generated by: LCOV version 2.0-1