Rationalized and tidied up polynomial functions.

No significant implications for program functionality.

SConstruct

@@ -69,7 +69,7 @@ variants['build_corfun']               = (['corfun.f90','maths_module.f90'],env_
 variants['build_diffusion']            = (['diffusion.f90','config_io_module.f90'],env_normal)
 variants['build_diffusion_test']       = (['diffusion_test.f90']+utnoavrgs,env_normal)
 variants['build_eos_lj']               = (['eos_lj.f90','eos_lj_module.f90','lrc_lj_module.f90'],env_normal)
-variants['build_eos_hs']               = (['eos_hs.f90'],env_normal)
+variants['build_eos_hs']               = (['eos_hs.f90','maths_module.f90'],env_normal)
 variants['build_error_calc']           = (['error_calc.f90','maths_module.f90'],env_normal)
 variants['build_ewald']                = (['ewald.f90','ewald_module.f90','mesh_module.f90','config_io_module.f90'],env_fftw)
 variants['build_fft3dwrap']            = (['fft3dwrap.f90'],env_fftw)

bd_nvt_lj.f90

@@ -187,6 +187,7 @@ SUBROUTINE b_propagator ( t )
   END SUBROUTINE b_propagator
 
   SUBROUTINE o_propagator ( t )
+    USE maths_module, ONLY : expm1
     IMPLICIT NONE
     REAL, INTENT(in) :: t ! Time over which to propagate (typically dt)
 
@@ -195,14 +196,8 @@ SUBROUTINE o_propagator ( t )
     REAL                 :: x, c
     REAL, DIMENSION(3,n) :: zeta
 
-    REAL, PARAMETER :: c1 = 2.0, c2 = -2.0, c3 = 4.0/3.0, c4 = -2.0/3.0 ! Taylor series coefficients
-
     x = gamma * t
-    IF ( ABS(x) > 0.0001 ) THEN ! Use formula
-       c = 1-EXP(-2.0*x)
-    ELSE ! Use Taylor expansion for low x
-       c = x * ( c1 + x * ( c2 + x * ( c3 + x * c4 ) ) ) 
-    END IF
+    c = -expm1(-2*x) ! 1-exp(-2*x)
     c = SQRT ( c )
 
     CALL random_normals ( 0.0, SQRT(temperature), zeta ) ! Random momenta

corfun.f90

@@ -41,7 +41,7 @@ PROGRAM corfun
 
   USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
   USE, INTRINSIC :: iso_c_binding
-  USE               maths_module,    ONLY : random_normal
+  USE               maths_module,    ONLY : random_normal, expm1
 
   IMPLICIT NONE
   INCLUDE 'fftw3.f03'
@@ -78,10 +78,6 @@ PROGRAM corfun
   COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE :: fft_out  ! data to be transformed (0:fft_len-1)
   TYPE(C_PTR)                                          :: fft_plan ! plan needed for FFTW
 
-  ! Taylor series coefficients
-  REAL, PARAMETER :: b1 = 2.0, b2 = -2.0,     b3 = 4.0/3.0, b4 = -2.0/3.0  ! b = 1.0 - EXP(-2.0*x)
-  REAL, PARAMETER :: d1 = 1.0, d2 = -1.0/2.0, d3 = 1.0/6.0, d4 = -1.0/24.0 ! d = 1.0 - EXP(-x)
-
   NAMELIST /nml/ nt, origin_interval, nstep, nequil, delta, temperature
 
   WRITE( unit=output_unit, fmt='(a)' ) 'corfun'
@@ -136,15 +132,10 @@ PROGRAM corfun
   ! Coefficients used in algorithm
   x = delta*kappa
   e = EXP(-x) ! theta in B&B paper
-  IF ( ABS(x) > 0.0001 ) THEN
-     b = 1.0 - EXP(-2.0*x)
-     d = 1.0 - EXP(-x)
-  ELSE  ! Taylor expansions for low x
-     b = x * ( b1 + x * ( b2 + x * ( b3 + x * b4 ) ) )
-     d = x * ( d1 + x * ( d2 + x * ( d3 + x * d4 ) ) )
-  END IF
-  b      = SQRT ( b )
-  b      = b * SQRT ( kappa/2.0 ) ! alpha in B&B paper
+  b = -expm1(-2*x) ! 1-exp(-2*x)
+  d = -expm1(-x)   ! 1-exp(-x)
+  b = SQRT ( b )
+  b = b * SQRT ( kappa/2.0 ) ! alpha in B&B paper
   stddev = SQRT(2.0*temperature)
 
   ! Data generation

diffusion_test.f90

@@ -139,18 +139,14 @@ SUBROUTINE a_propagator ( t ) ! A propagator (drift)
   END SUBROUTINE a_propagator
 
   SUBROUTINE o_propagator ( t ) ! O propagator (friction and random contributions)
+    USE maths_module, ONLY : expm1
     IMPLICIT NONE
     REAL, INTENT(in) :: t ! Time over which to propagate (typically dt)
 
     REAL            :: x, c
-    REAL, PARAMETER :: c1 = 2.0, c2 = -2.0, c3 = 4.0/3.0, c4 = -2.0/3.0
 
     x = gamma * t
-    IF ( ABS(x) > 0.0001 ) THEN
-       c = 1-EXP(-2.0*x)
-    ELSE
-       c = x * ( c1 + x * ( c2 + x * ( c3 + x * c4 )) ) ! Taylor expansion for low x
-    END IF
+    c = -expm1(-2*x) ! 1-exp(-2*x)
     c = SQRT ( c )
 
     CALL random_normals ( 0.0, SQRT(temperature), zeta ) ! Random momenta

dpd_module.f90

@@ -314,6 +314,7 @@ SUBROUTINE shardlow ( box, temperature, gamma_step )
   END SUBROUTINE shardlow
 
   FUNCTION p_approx ( a, rho, temperature ) RESULT ( p )
+    USE maths_module, ONLY : polyval
     REAL             :: p           ! Returns approximate pressure
     REAL, INTENT(in) :: a           ! Force strength parameter
     REAL, INTENT(in) :: rho         ! Density
@@ -324,30 +325,17 @@ FUNCTION p_approx ( a, rho, temperature ) RESULT ( p )
     REAL, DIMENSION(0:4), PARAMETER :: b = [ 9.755e-2, -4.912e-3, 1.556e-4, -2.585e-6, 1.705e-8 ]
     REAL,                 PARAMETER :: c1 = 0.0802, c2 = 0.7787
 
-    ! This expression is given by Groot and Warren, J Chem Phys 107, 4423 (1997)
-!    alpha = 0.101
+    ! Original expression given by Groot and Warren, J Chem Phys 107, 4423 (1997), alpha = 0.101
 
     ! This is the revised formula due to Liyana-Arachchi, Jamadagni, Elke, Koenig, Siepmann,
     ! J Chem Phys 142, 044902 (2015)
-    b2    = polynomial ( b, a )                                         ! This is B2/a, eqn (10) of above paper
-    alpha = b2 / ( 1.0 + rho**3 ) + ( c1*rho**2 ) / ( 1.0 + c2*rho**2 ) ! This is eqn (14) of above paper
+    ! Eqn (10) of above paper, B2/a
+    b2 = polyval ( a, b )
+    ! Eqn (14) of above paper
+    alpha = b2 / ( 1.0 + rho**3 ) + ( c1*rho**2 ) / ( 1.0 + c2*rho**2 )
 
     p = rho * temperature + alpha * a * rho**2
 
   END FUNCTION p_approx
 
-  FUNCTION polynomial ( c, x ) RESULT ( f )
-    REAL, DIMENSION(:), INTENT(in) :: c ! coefficients of x**0, x**1, x**2 etc
-    REAL,               INTENT(in) :: x ! argument
-    REAL                           :: f ! result
-
-    INTEGER :: i
-
-    f = 0.0
-    DO i = SIZE(c), 1, -1
-       f = f * x + c(i)
-    END DO
-
-  END FUNCTION polynomial
-
 END MODULE dpd_module

eos_hs.f90

@@ -32,6 +32,8 @@ PROGRAM eos_hs
 
   USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
 
+  USE maths_module, ONLY : polyval
+
   IMPLICIT NONE
 
   REAL            :: density, p, z, eta
@@ -62,25 +64,9 @@ PROGRAM eos_hs
 
   ! Equation (6) of Hansen-Goos (2016)
   z = a * LOG ( 1.0-eta ) / eta
-  z = z + polynomial ( b, eta ) / ( 1.0 - eta ) ** 3 ! Compressibility factor P/(rho*kT)
-  p = z * density                                    ! Pressure P / kT
+  z = z + polyval ( eta, b ) / ( 1.0 - eta ) ** 3 ! Compressibility factor P/(rho*kT)
+  p = z * density                                 ! Pressure P / kT
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Pressure P',                            p
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Compressibility factor Z = P/(rho*kT)', z
 
-CONTAINS
-
-  FUNCTION polynomial ( c, x ) RESULT ( f )
-    REAL, DIMENSION(:), INTENT(in) :: c ! coefficients of x**0, x**1, x**2 etc
-    REAL,               INTENT(in) :: x ! argument
-    REAL                           :: f ! result
-
-    INTEGER :: i
-
-    f = 0.0
-    DO i = SIZE(c), 1, -1
-       f = f * x + c(i)
-    END DO
-
-  END FUNCTION polynomial
-
 END PROGRAM eos_hs

error_calc.f90

@@ -32,7 +32,7 @@ PROGRAM error_calc
   ! and AD Baczewski and SD Bond J Chem Phys 139 044107 (2013)
 
   USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
-  USE               maths_module,    ONLY : random_normal
+  USE               maths_module,    ONLY : random_normal, expm1
 
   IMPLICIT NONE
 
@@ -51,10 +51,6 @@ PROGRAM error_calc
   REAL    :: average, variance, stddev, err, si, tcor, x, e, b, d
   REAL    :: a_blk, a_run, a_avg, a_var, a_var_1, a_err
 
-  ! Taylor series coefficients
-  REAL, PARAMETER :: b1 = 2.0, b2 = -2.0,     b3 = 4.0/3.0, b4 = -2.0/3.0  ! b = 1.0 - EXP(-2.0*x)
-  REAL, PARAMETER :: d1 = 1.0, d2 = -1.0/2.0, d3 = 1.0/6.0, d4 = -1.0/24.0 ! d = 1.0 - EXP(-x)
-
   NAMELIST /nml/ nstep, nequil, n_repeat, delta, variance, average
 
   ! Example default values
@@ -97,15 +93,10 @@ PROGRAM error_calc
   ! Coefficients used in algorithm
   x = delta*kappa
   e = EXP(-x) ! theta in B&B paper
-  IF ( ABS(x) > 0.0001 ) THEN
-     b = 1.0 - EXP(-2.0*x)
-     d = 1.0 - EXP(-x)
-  ELSE ! Taylor expansions for low x
-     b = x * ( b1 + x * ( b2 + x * ( b3 + x * b4 ) ) ) 
-     d = x * ( d1 + x * ( d2 + x * ( d3 + x * d4 ) ) )
-  END IF
-  b      = SQRT ( b )
-  b      = b * SQRT ( kappa/2.0 ) ! alpha in B&B paper  
+  b = -expm1(-2*x) ! 1-exp(-2*x)
+  d = -expm1(-x)   ! 1-exp(-x)
+  b = SQRT ( b )
+  b = b * SQRT ( kappa/2.0 ) ! alpha in B&B paper  
   stddev = SQRT(2.0*variance)     ! NB stddev of random forces, not data
 
   ! For this process, the results of interest can be calculated exactly

maths_module.f90

@@ -41,7 +41,8 @@ MODULE maths_module
   PUBLIC :: metropolis
 
   ! Public low-level mathematical routines and string operations
-  PUBLIC :: rotate_vector, rotate_quaternion, cross_product, outer_product, q_to_a, lowercase, solve
+  PUBLIC :: rotate_vector, rotate_quaternion, cross_product, outer_product, q_to_a
+  PUBLIC :: lowercase, solve, polyval, expm1, expm1o
 
   ! Public order parameter calculations
   PUBLIC :: orientational_order, translational_order, nematic_order
@@ -902,4 +903,64 @@ FUNCTION q_to_a ( q ) RESULT ( a )
     a(3,:) = [     2*(q(1)*q(3)+q(0)*q(2)),       2*(q(2)*q(3)-q(0)*q(1)),   q(0)**2-q(1)**2-q(2)**2+q(3)**2 ] ! 3rd row
   END FUNCTION q_to_a
 
+  FUNCTION polyval ( x, c ) RESULT ( f )
+    IMPLICIT NONE
+    REAL                            :: f ! Returns polynomial in ...
+    REAL,                INTENT(in) :: x ! argument
+    REAL, DIMENSION(0:), INTENT(in) :: c ! given coefficients (ascending powers of x)
+
+    ! Uses Horner's rule
+
+    INTEGER :: i
+
+    f = c(UBOUND(c,1))
+    DO i = UBOUND(c,1)-1, 0, -1
+       f = f * x + c(i)
+    END DO
+  END FUNCTION polyval
+
+  FUNCTION expm1 ( x ) RESULT ( f )
+    IMPLICIT NONE
+    REAL, INTENT(in) :: x ! Argument
+    REAL             :: f ! Returns value of exp(x)-1
+
+    ! At small x, imprecision may arise through subtraction of two values O(1).
+    ! There are various ways of avoiding this.
+    ! We follow some others and use the identity: exp(x)-1 = x*exp(x/2)*[sinh(x/2)/(x/2)].
+    ! For small x, sinh(x)/x = g0 + g1*x**2 + g2*x**4 + ...
+    ! where the coefficient of x**(2n) is gn = 1/(2*n+1)!
+
+    REAL, DIMENSION(0:4), PARAMETER :: g = 1.0 / [1,6,120,5040,362880]
+    REAL,                 PARAMETER :: tol = 0.01
+
+    IF ( ABS(x) > tol ) THEN
+       f = EXP(x) - 1.0
+    ELSE
+       f = x * EXP(x/2) * polyval ( (x/2)**2, g )
+    END IF
+
+  END FUNCTION expm1
+
+  FUNCTION expm1o ( x ) RESULT ( f )
+    IMPLICIT NONE
+    REAL, INTENT(in) :: x ! Argument
+    REAL             :: f ! Returns value of (exp(x)-1)/x
+
+    ! At small x, we must guard against the ratio of imprecise small values.
+    ! There are various ways of avoiding this.
+    ! We follow some others and use the identity: (exp(x)-1)/x = exp(x/2)*[sinh(x/2)/(x/2)].
+    ! For small x, sinh(x)/x = g0 + g1*x**2 + g2*x**4 + ...
+    ! where the coefficient of x**(2n) is gn = 1/(2*n+1)!
+
+    REAL, DIMENSION(0:4), PARAMETER :: g = 1.0 / [1,6,120,5040,362880]
+    REAL,                 PARAMETER :: tol = 0.01
+
+    IF ( ABS(x) > tol ) THEN
+       f = ( EXP(x) - 1.0 ) / x
+    ELSE
+       f = EXP(x/2) * polyval ( (x/2)**2, g )
+    END IF
+
+  END FUNCTION expm1o
+
 END MODULE maths_module

md_npt_lj.f90

@@ -222,24 +222,19 @@ PROGRAM md_npt_lj
 CONTAINS
 
   SUBROUTINE u1_propagator ( t ) ! U1 and U1' combined: position and strain drift propagator
+    USE maths_module, ONLY : expm1o
     IMPLICIT NONE
     REAL, INTENT(in) :: t ! Time over which to propagate (typically dt)
 
-    REAL            :: x, c
-    REAL, PARAMETER :: c1 = -1.0/2.0, c2 = 1.0/6.0, c3 = -1.0/24.0
+    REAL :: x, c
 
     ! U1 part
     ! The propagator for r(:,:) looks different from the formula given on p142 of the text.
     ! However it is easily derived from that formula, bearing in mind that coordinates in
     ! this program are divided by the box length, which is itself updated in this routine.
 
     x = t * p_eps / w_eps ! Time step * time derivative of strain
-
-    IF ( ABS(x) < 0.001 ) THEN ! Guard against small values
-       c = 1.0 + x * ( c1 + x * ( c2 + x * c3 ) ) ! Taylor series to order 3
-    ELSE
-       c = (1.0-EXP(-x))/x
-    END IF ! End guard against small values
+    c = expm1o(-x) ! (1-exp(-x))/x
 
     r(:,:) = r(:,:) + c * t * v(:,:) / box ! Positions in box=1 units
     r(:,:) = r(:,:) - ANINT ( r(:,:) )     ! Periodic boundaries
@@ -254,20 +249,15 @@ SUBROUTINE u1_propagator ( t ) ! U1 and U1' combined: position and strain drift
   END SUBROUTINE u1_propagator
 
   SUBROUTINE u2_propagator ( t ) ! U2: velocity kick step propagator
+    USE maths_module, ONLY : expm1o
     IMPLICIT NONE
     REAL, INTENT(in) :: t ! Time over which to propagate (typically dt/2)
 
-    REAL            :: x, c, alpha
-    REAL, PARAMETER :: c1 = -1.0/2.0, c2 = 1.0/6.0, c3 = -1.0/24.0
+    REAL :: x, c, alpha
 
     alpha = 1.0 + 3.0 / g
     x = t * alpha * p_eps / w_eps
-
-    IF ( ABS(x) < 0.001 ) THEN ! Guard against small values
-       c = 1.0 + x * ( c1 + x * ( c2 + x * c3 ) ) ! Taylor series to order 3
-    ELSE
-       c = (1.0-EXP(-x))/x
-    END IF ! End guard against small values
+    c = expm1o(-x) ! (1-exp(-x))/x
 
     v(:,:) = v(:,:)*EXP(-x) + c * t * f(:,:)
 
@@ -302,13 +292,13 @@ SUBROUTINE u3_propagator ( t ) ! U3 and U3' combined: thermostat propagator
   END SUBROUTINE u3_propagator
 
   SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4 and U4' combined: thermostat propagator
+    USE maths_module, ONLY : expm1o
     IMPLICIT NONE
     REAL,    INTENT(in) :: t               ! Time over which to propagate (typically dt/4)
     INTEGER, INTENT(in) :: j_start, j_stop ! Order in which to tackle variables
 
-    INTEGER         :: j, j_stride
-    REAL            :: gj, x, c
-    REAL, PARAMETER :: c1 = -1.0/2.0, c2 = 1.0/6.0, c3 = -1.0/24.0
+    INTEGER :: j, j_stride
+    REAL    :: gj, x, c
 
     IF ( j_start > j_stop ) THEN
        j_stride = -1
@@ -333,12 +323,7 @@ SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4 and U4' combined: thermosta
        ELSE
 
           x = t * p_eta(j+1)/q(j+1)
-
-          IF ( ABS(x) < 0.001 ) THEN ! Guard against small values
-             c = 1.0 + x * ( c1 + x * ( c2 + x * c3 ) ) ! Taylor series to order 3
-          ELSE
-             c = (1.0-EXP(-x))/x
-          END IF ! End guard against small values
+          c = expm1o(-x) ! (1-exp(-x))/x
 
           p_eta(j) = p_eta(j)*EXP(-x) + t * gj * c
 
@@ -363,12 +348,7 @@ SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4 and U4' combined: thermosta
        ELSE
 
           x = t * p_eta_baro(j+1)/q_baro(j+1)
-
-          IF ( ABS(x) < 0.001 ) THEN ! Guard against small values
-             c = 1.0 + x * ( c1 + x * ( c2 + x * c3 ) ) ! Taylor series to order 3
-          ELSE
-             c = (1.0-EXP(-x))/x
-          END IF ! End guard against small values
+          c = expm1o(-x) ! (1-exp(-x))/x
 
           p_eta_baro(j) = p_eta_baro(j)*EXP(-x) + t * gj * c