From 6e5eea87b75ac69bba8c76946d652c0f2efd47c1 Mon Sep 17 00:00:00 2001 From: Michael-P-Allen Date: Sun, 21 Jan 2024 17:31:01 +0000 Subject: [PATCH] Rationalized and tidied up polynomial functions. No significant implications for program functionality. --- SConstruct | 2 +- bd_nvt_lj.f90 | 9 +---- corfun.f90 | 19 +++------- diffusion_test.f90 | 8 +--- dpd_module.f90 | 24 +++--------- eos_hs.f90 | 22 ++--------- error_calc.f90 | 19 +++------- maths_module.f90 | 63 ++++++++++++++++++++++++++++++- md_npt_lj.f90 | 42 ++++++--------------- md_nvt_lj.f90 | 13 ++----- python_examples/bd_nvt_lj.py | 2 +- python_examples/corfun.py | 13 ++----- python_examples/diffusion_test.py | 6 +-- python_examples/dpd_module.py | 15 +++++--- python_examples/eos_hs.py | 7 ++-- python_examples/error_calc.py | 13 ++----- python_examples/maths_module.py | 18 +++++++++ python_examples/md_npt_lj.py | 11 ++++-- python_examples/md_nvt_lj.py | 3 +- python_examples/qmc_pi_sho.py | 27 ++++++------- qmc_pi_sho.f90 | 40 +++++++------------- 21 files changed, 181 insertions(+), 195 deletions(-) diff --git a/SConstruct b/SConstruct index 90037b8..2abc39b 100644 --- a/SConstruct +++ b/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) diff --git a/bd_nvt_lj.f90 b/bd_nvt_lj.f90 index ca79d7c..bc7805f 100644 --- a/bd_nvt_lj.f90 +++ b/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 diff --git a/corfun.f90 b/corfun.f90 index cf93c19..5b6103d 100644 --- a/corfun.f90 +++ b/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 diff --git a/diffusion_test.f90 b/diffusion_test.f90 index 28def4b..fab91ed 100644 --- a/diffusion_test.f90 +++ b/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 diff --git a/dpd_module.f90 b/dpd_module.f90 index 6b4285e..a2f9e21 100644 --- a/dpd_module.f90 +++ b/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 diff --git a/eos_hs.f90 b/eos_hs.f90 index f627692..cbe087e 100644 --- a/eos_hs.f90 +++ b/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 diff --git a/error_calc.f90 b/error_calc.f90 index 84504d1..da5f40e 100644 --- a/error_calc.f90 +++ b/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 diff --git a/maths_module.f90 b/maths_module.f90 index 57c46d3..7a430b1 100644 --- a/maths_module.f90 +++ b/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 diff --git a/md_npt_lj.f90 b/md_npt_lj.f90 index a691a22..3d73707 100644 --- a/md_npt_lj.f90 +++ b/md_npt_lj.f90 @@ -222,11 +222,11 @@ 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. @@ -234,12 +234,7 @@ SUBROUTINE u1_propagator ( t ) ! U1 and U1' combined: position and strain drift ! 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 diff --git a/md_nvt_lj.f90 b/md_nvt_lj.f90 index f0cc77c..8e0d0c4 100644 --- a/md_nvt_lj.f90 +++ b/md_nvt_lj.f90 @@ -221,13 +221,13 @@ SUBROUTINE u3_propagator ( t ) ! U3: thermostat propagator END SUBROUTINE u3_propagator SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4: 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 @@ -250,12 +250,7 @@ SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4: thermostat propagator 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 diff --git a/python_examples/bd_nvt_lj.py b/python_examples/bd_nvt_lj.py index 56bc62a..893d677 100755 --- a/python_examples/bd_nvt_lj.py +++ b/python_examples/bd_nvt_lj.py @@ -122,7 +122,7 @@ def o_propagator ( t ): import numpy as np x = gamma*t - c = 1-np.exp(-2*x) if np.fabs(x) > 0.0001 else np.polyval([-2/3,4/3,-2.0,2.0,0.0],x) + c = -np.expm1(-2*x) # 1-exp(-2*x) c = np.sqrt(c) v = v*np.exp(-x) + c*np.sqrt(temperature)*np.random.randn(n,3) diff --git a/python_examples/corfun.py b/python_examples/corfun.py index 7317762..b78567e 100755 --- a/python_examples/corfun.py +++ b/python_examples/corfun.py @@ -113,15 +113,10 @@ def c_exact(t): # Coefficients used in algorithm x = delta*kappa e = math.exp(-x) # theta in B&B paper -if math.fabs(x) > 0.0001: - b = 1.0 - math.exp(-2.0*x) - d = 1.0 - math.exp(-x) -else: # Taylor expansions for low x (coefficients of decreasing powers) - b = np.polyval([-2.0/3.0,4.0/3.0,-2.0,2.0,0.0],x) - d = np.polyval([-1.0/24.0,1.0/6.0,-1.0/2.0,1.0,0.0],x) - -b = math.sqrt ( b ) -b = b * math.sqrt ( kappa/2.0 ) # alpha in B&B paper +b = -np.expm1(-2*x) # 1-exp(-2*x) +d = -np.expm1(-x) # 1-exp(-x) +b = math.sqrt ( b ) +b = b * math.sqrt ( kappa/2.0 ) # alpha in B&B paper stddev = math.sqrt(2.0*temperature) # Data generation diff --git a/python_examples/diffusion_test.py b/python_examples/diffusion_test.py index 2f11fff..22ca8a9 100755 --- a/python_examples/diffusion_test.py +++ b/python_examples/diffusion_test.py @@ -43,13 +43,11 @@ def o_propagator ( t, v ): gamma and temperature are accessed from the calling program The function returns the new velocities. """ - import math import numpy as np x = gamma*t - # Taylor expansion for small x - c = 1.0-math.exp(-2.0*x) if math.fabs(x)>0.0001 else np.polyval([-2.0/3.0,4.0/3.0,-2.0,2.0,0.0],x) - c = math.sqrt(c) + c = -np.expm1(-2*x) # 1-exp(-2*x) + c = np.sqrt(c) return np.exp(-x) * v + c * np.sqrt(temperature) * np.random.randn(n,3) # diffusion_test program diff --git a/python_examples/dpd_module.py b/python_examples/dpd_module.py index 7eb7765..ff4b6e7 100644 --- a/python_examples/dpd_module.py +++ b/python_examples/dpd_module.py @@ -195,14 +195,19 @@ def shardlow ( box, temperature, gamma_step, v, pairs ): def p_approx ( a, rho, temperature ): """Returns approximate pressure.""" + import numpy as np + from numpy.polynomial.polynomial import polyval + # 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) - import numpy as np - c1, c2 = 0.0802, 0.7787 - b = [1.705e-8,-2.585e-6,1.556e-4,-4.912e-3,9.755e-2] - b2 = np.polyval (b,a) # This is B2/a, eqn (10) of above paper - alpha = b2 / ( 1 + rho**3 ) + ( c1*rho**2 ) / ( 1 + c2*rho**2 ) # This is eqn (14) of above paper + b = [ 9.755e-2, -4.912e-3, 1.556e-4, -2.585e-6, 1.705e-8 ] + # Eqn (10) of above paper, B2/a + b2 = polyval (a,b) + # Eqn (14) of above paper + alpha = b2 / ( 1 + rho**3 ) + ( c1*rho**2 ) / ( 1 + c2*rho**2 ) + return rho * temperature + alpha * a * rho**2 diff --git a/python_examples/eos_hs.py b/python_examples/eos_hs.py index 5a81eaf..f39bf57 100755 --- a/python_examples/eos_hs.py +++ b/python_examples/eos_hs.py @@ -29,6 +29,7 @@ import json import sys import numpy as np +from numpy.polynomial.polynomial import polyval # This program uses the function derived by H Hansen-Goos, J Chem Phys, 16, 164506 (2016) # which is claimed to be an excellent fit to simulation data over the whole fluid density range @@ -38,7 +39,7 @@ # The coefficients appear in Table I of Hansen-Goos (2016) a = 8.0 -b = [ -0.096495, 0.248245, 0.041212, -1.265575, -2.635232, 47.0/3.0, -19.0, 9.0 ] +b = [ 9.0, -19.0, 47.0/3.0, -2.635232, -1.265575, 0.041212, 0.248245, -0.096495 ] # Read parameters in JSON format try: @@ -65,7 +66,7 @@ # Equation (6) of Hansen-Goos (2016) z = a * np.log ( 1.0-eta ) / eta -z = z + np.polyval ( 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 print ( "{:40}{:15.6f}".format('Pressure P', p ) ) print ( "{:40}{:15.6f}".format('Compressibility factor Z = P/(rho*kT)', z ) ) diff --git a/python_examples/error_calc.py b/python_examples/error_calc.py index 8814f8b..8d09b81 100755 --- a/python_examples/error_calc.py +++ b/python_examples/error_calc.py @@ -84,15 +84,10 @@ # Coefficients used in algorithm x = delta*kappa e = math.exp(-x) # theta in B&B paper -if math.fabs(x) > 0.0001: - b = 1.0 - math.exp(-2.0*x) - d = 1.0 - math.exp(-x) -else: # Taylor expansions for low x (coefficients of decreasing powers) - b = np.polyval([-2.0/3.0,4.0/3.0,-2.0,2.0,0.0],x) - d = np.polyval([-1.0/24.0,1.0/6.0,-1.0/2.0,1.0,0.0],x) - -b = math.sqrt ( b ) -b = b * math.sqrt ( kappa/2.0 ) # alpha in B&B paper +b = -np.expm1(-2*x) # 1-exp(-2*x) +d = -np.expm1(-x) # 1-exp(-x) +b = math.sqrt ( b ) +b = b * math.sqrt ( kappa/2.0 ) # alpha in B&B paper stddev = math.sqrt(2.0*variance) # NB stddev of random forces, not data # For this process, the results of interest can be calculated exactly diff --git a/python_examples/maths_module.py b/python_examples/maths_module.py index b88b52d..ad5c6f7 100644 --- a/python_examples/maths_module.py +++ b/python_examples/maths_module.py @@ -249,3 +249,21 @@ def q_to_a ( q ): return a +def expm1o ( x ): + """Returns (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)! + # NB for the related function exp(x)-1 we rely on the built-in expm1(x) + + import numpy as np + from numpy.polynomial.polynomial import polyval + g = 1.0 / np.array([1,6,120,5040,362880],dtype=np.float_) + tol = 0.01 + if abs(x) > tol: + return np.expm1(x)/x + else: + return np.exp(x/2)*polyval((x/2)**2,g) diff --git a/python_examples/md_npt_lj.py b/python_examples/md_npt_lj.py index fb9c716..8bcdef6 100755 --- a/python_examples/md_npt_lj.py +++ b/python_examples/md_npt_lj.py @@ -106,6 +106,7 @@ def u1_propagator ( t ): global r, eps, box import numpy as np + from maths_module import expm1o # U1 part # The propagator for r looks different from the formula given on p142 of the text. @@ -113,7 +114,7 @@ def u1_propagator ( t ): # 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 - c = (1.0-np.exp(-x))/x if np.fabs(x)>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values + c = expm1o(-x) # (1-exp(-x))/x r = r + c * t * v / box # Positions in box=1 units r = r - np.rint ( r ) # Periodic boundaries @@ -134,10 +135,11 @@ def u2_propagator ( t ): global v import numpy as np + from maths_module import expm1o alpha = 1.0 + 3 / g x = t * alpha * p_eps / w_eps - c = (1.0-np.exp(-x))/x if np.fabs(x)>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values + c = expm1o(-x) # (1-exp(-x))/x v = v*np.exp(-x) + c * t * f @@ -185,6 +187,7 @@ def u4_propagator ( t, j_list ): global p_eta, p_eta_baro import numpy as np + from maths_module import expm1o # U4 part @@ -197,7 +200,7 @@ def u4_propagator ( t, j_list ): p_eta[j] = p_eta[j] + t * gj # The equation for p_eta[M-1] is different else: x = t * p_eta[j+1]/q[j+1] - c = (1.0-np.exp(-x))/x if np.fabs(x)>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values + c = expm1o(-x) # (1-exp(-x))/x p_eta[j] = p_eta[j]*np.exp(-x) + t * gj * c # U4' part @@ -211,7 +214,7 @@ def u4_propagator ( t, j_list ): p_eta_baro[j] = p_eta_baro[j] + t * gj # The equation for p_eta_baro[M-1] is different else: x = t * p_eta_baro[j+1]/q_baro[j+1] - c = (1.0-np.exp(-x))/x if np.fabs(x)>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values + c = expm1o(-x) # (1-exp(-x))/x p_eta_baro[j] = p_eta_baro[j]*np.exp(-x) + t * gj * c # Takes in a configuration of atoms (positions, velocities) diff --git a/python_examples/md_nvt_lj.py b/python_examples/md_nvt_lj.py index 165f80b..e92a8fc 100755 --- a/python_examples/md_nvt_lj.py +++ b/python_examples/md_nvt_lj.py @@ -144,6 +144,7 @@ def u4_propagator ( t, j_list ): global p_eta import numpy as np + from maths_module import expm1o for j in j_list: if j==0: @@ -154,7 +155,7 @@ def u4_propagator ( t, j_list ): p_eta[j] = p_eta[j] + t * gj # The equation for p_eta[M-1] is different else: x = t * p_eta[j+1]/q[j+1] - c = (1.0-np.exp(-x))/x if np.fabs(x)>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values + c = expm1o(-x) # (1-exp(-x))/x p_eta[j] = p_eta[j]*np.exp(-x) + t * gj * c # Takes in a configuration of atoms (positions, velocities) diff --git a/python_examples/qmc_pi_sho.py b/python_examples/qmc_pi_sho.py index d62b291..8fc00fd 100755 --- a/python_examples/qmc_pi_sho.py +++ b/python_examples/qmc_pi_sho.py @@ -68,7 +68,8 @@ def e_pi_sho ( p, beta ): import math import numpy as np - + from numpy.polynomial.polynomial import polyval + assert p>0, 'Error in value of p' t = 1 / beta @@ -77,30 +78,30 @@ def e_pi_sho ( p, beta ): if p==1: e = t elif p==2: - e = 1.0 + 1.0 / np.polyval([4,1],s) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) e = e * t elif p==3: - e = 1.0 + 2.0 / np.polyval([3,1],s) + e = 1.0 + 2.0 / polyval ( s, [1.,3.] ) e = e * t elif p==4: - e = 1.0 + 1.0 / np.polyval([4,1],s) - e = e + 2.0 / np.polyval([2,1],s) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) + e = e + 2.0 / polyval ( s, [1.,2.] ) e = e * t elif p==5: - e = 1.0 + np.polyval([10,4],s)/np.polyval([5,5,1],s) + e = 1.0 + polyval ( s, [4.,10.] )/polyval ( s, [1.,5.,5.] ) e = e * t elif p==6: - e = 1.0 + 1.0 / np.polyval([4,1],s) - e = e + 2.0 / np.polyval([1,1],s) - e = e + 2.0 / np.polyval([3,1],s) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) + e = e + 2.0 / polyval ( s, [1.,1.] ) + e = e + 2.0 / polyval ( s, [1.,3.] ) e = e * t elif p==7: - e = 1.0 + np.polyval([28,28,6],s) / np.polyval([7,14,7,1],s) + e = 1.0 + polyval ( s, [6.,28.,28.] ) / polyval ( s, [1.,7.,14.,7.] ) e = e * t elif p==8: - e = 1.0 + 1.0 / np.polyval([4,1],s) - e = e + 2.0 / np.polyval([2,1],s) - e = e + np.polyval ([8,4],s) / np.polyval([2,4,1],s) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) + e = e + 2.0 / polyval ( s, [1.,2.] ) + e = e + polyval ( s, [4.,8.] ) / polyval ( s, [1.,4.,2.] ) e = e * t else: alpha = 0.5 * beta / p diff --git a/qmc_pi_sho.f90 b/qmc_pi_sho.f90 index 5b93527..830d584 100644 --- a/qmc_pi_sho.f90 +++ b/qmc_pi_sho.f90 @@ -202,6 +202,7 @@ FUNCTION calc_variables ( ) RESULT ( variables ) END FUNCTION calc_variables FUNCTION e_pi_sho ( p, beta ) RESULT ( e ) + USE maths_module, ONLY : polyval IMPLICIT NONE INTEGER, INTENT(in) :: p REAL, INTENT(in) :: beta @@ -234,36 +235,36 @@ FUNCTION e_pi_sho ( p, beta ) RESULT ( e ) e = t CASE (2) - e = 1.0 + 1.0 / polynomial ( [1,4], s ) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) e = e * t CASE (3) - e = 1.0 + 2.0 / polynomial ( [1,3], s ) + e = 1.0 + 2.0 / polyval ( s, [1.,3.] ) e = e * t CASE (4) - e = 1.0 + 1.0 / polynomial ( [1,4], s ) - e = e + 2.0 / polynomial ( [1,2], s ) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) + e = e + 2.0 / polyval ( s, [1.,2.] ) e = e * t CASE (5) - e = 1.0 + polynomial ( [4,10], s ) / polynomial ( [1,5,5], s ) + e = 1.0 + polyval ( s, [4.,10.] ) / polyval ( s, [1.,5.,5.] ) e = e * t CASE (6) - e = 1.0 + 1.0 / polynomial ( [1,4], s ) - e = e + 2.0 / polynomial ( [1,1], s ) - e = e + 2.0 / polynomial ( [1,3], s ) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) + e = e + 2.0 / polyval ( s, [1.,1.] ) + e = e + 2.0 / polyval ( s, [1.,3.] ) e = e * t CASE (7) - e = 1.0 + polynomial ( [6,28,28], s ) / polynomial ( [1,7,14,7], s ) + e = 1.0 + polyval ( s, [6.,28.,28.] ) / polyval ( s, [1.,7.,14.,7.] ) e = e * t CASE (8) - e = 1.0 + 1.0 / polynomial ( [1,4], s ) - e = e + 2.0 / polynomial ( [1,2], s ) - e = e + polynomial ( [4,8], s ) / polynomial ( [1,4,2], s ) + e = 1.0 + 1.0 / polyval ( s, [1.,4.] ) + e = e + 2.0 / polyval ( s, [1.,2.] ) + e = e + polyval ( s, [4.,8.] ) / polyval ( s, [1.,4.,2.] ) e = e * t CASE default @@ -277,20 +278,5 @@ FUNCTION e_pi_sho ( p, beta ) RESULT ( e ) END FUNCTION e_pi_sho - FUNCTION polynomial ( c, x ) RESULT ( f ) - IMPLICIT NONE - INTEGER, 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 + REAL ( c(i) ) - END DO - - END FUNCTION polynomial - END PROGRAM qmc_pi_sho