Corrected x -> abs(x) to select various low-x polynomial approximants.

In principle this only affects the md_nvt_lj and md_npt_lj codes,
and in practice the effect is imperceptible since the approximants
are very accurate in the applicable range of x.

bd_nvt_lj.f90

@@ -198,7 +198,7 @@ SUBROUTINE o_propagator ( t )
     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 ( x > 0.0001 ) THEN ! Use formula
+    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 ) ) ) 

corfun.f90

@@ -136,7 +136,7 @@ PROGRAM corfun
   ! Coefficients used in algorithm
   x = delta*kappa
   e = EXP(-x) ! theta in B&B paper
-  IF ( x > 0.0001 ) THEN
+  IF ( ABS(x) > 0.0001 ) THEN
      b = 1.0 - EXP(-2.0*x)
      d = 1.0 - EXP(-x)
   ELSE  ! Taylor expansions for low x

diffusion_test.f90

@@ -146,7 +146,7 @@ SUBROUTINE o_propagator ( t ) ! O propagator (friction and random contributions)
     REAL, PARAMETER :: c1 = 2.0, c2 = -2.0, c3 = 4.0/3.0, c4 = -2.0/3.0
 
     x = gamma * t
-    IF ( x > 0.0001 ) THEN
+    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

error_calc.f90

@@ -97,7 +97,7 @@ PROGRAM error_calc
   ! Coefficients used in algorithm
   x = delta*kappa
   e = EXP(-x) ! theta in B&B paper
-  IF ( x > 0.0001 ) THEN
+  IF ( ABS(x) > 0.0001 ) THEN
      b = 1.0 - EXP(-2.0*x)
      d = 1.0 - EXP(-x)
   ELSE ! Taylor expansions for low x

md_npt_lj.f90

@@ -235,7 +235,7 @@ SUBROUTINE u1_propagator ( t ) ! U1 and U1' combined: position and strain drift
 
     x = t * p_eps / w_eps ! Time step * time derivative of strain
 
-    IF ( x < 0.001 ) THEN ! Guard against small values
+    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
@@ -263,7 +263,7 @@ SUBROUTINE u2_propagator ( t ) ! U2: velocity kick step propagator
     alpha = 1.0 + 3.0 / g
     x = t * alpha * p_eps / w_eps
 
-    IF ( x < 0.001 ) THEN ! Guard against small values
+    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
@@ -334,7 +334,7 @@ SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4 and U4' combined: thermosta
 
           x = t * p_eta(j+1)/q(j+1)
 
-          IF ( x < 0.001 ) THEN ! Guard against small values
+          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
@@ -364,7 +364,7 @@ SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4 and U4' combined: thermosta
 
           x = t * p_eta_baro(j+1)/q_baro(j+1)
 
-          IF ( x < 0.001 ) THEN ! Guard against small values
+          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

md_nvt_lj.f90

@@ -251,7 +251,7 @@ SUBROUTINE u4_propagator ( t, j_start, j_stop ) ! U4: thermostat propagator
 
           x = t * p_eta(j+1)/q(j+1)
 
-          IF ( x < 0.001 ) THEN ! Guard against small values
+          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

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 x > 0.0001 else np.polyval([-2/3,4/3,-2.0,2.0,0.0],x)
+    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.sqrt(c)
     v = v*np.exp(-x) + c*np.sqrt(temperature)*np.random.randn(n,3)
 

python_examples/corfun.py

@@ -113,7 +113,7 @@ def c_exact(t):
 # Coefficients used in algorithm
 x = delta*kappa
 e = math.exp(-x) # theta in B&B paper
-if x > 0.0001:
+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)

python_examples/diffusion_test.py

@@ -48,7 +48,7 @@ def o_propagator ( t, v ):
 
     x = gamma*t
     # Taylor expansion for small x
-    c = 1.0-math.exp(-2.0*x) if x>0.0001 else np.polyval([-2.0/3.0,4.0/3.0,-2.0,2.0,0.0],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)
     return np.exp(-x) * v + c * np.sqrt(temperature) * np.random.randn(n,3)
 

python_examples/error_calc.py

@@ -84,7 +84,7 @@
 # Coefficients used in algorithm
 x = delta*kappa
 e = math.exp(-x) # theta in B&B paper
-if x > 0.0001:
+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)

python_examples/md_npt_lj.py

@@ -113,7 +113,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 x>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values
+    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
 
     r = r + c * t * v / box   # Positions in box=1 units
     r = r - np.rint ( r ) # Periodic boundaries
@@ -137,7 +137,7 @@ def u2_propagator ( t ):
 
     alpha = 1.0 + 3 / g
     x = t * alpha * p_eps / w_eps
-    c = (1.0-np.exp(-x))/x if x>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values
+    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
 
     v = v*np.exp(-x) + c * t * f
 
@@ -197,7 +197,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 x>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values
+            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
             p_eta[j] = p_eta[j]*np.exp(-x) + t * gj * c
 
     # U4' part
@@ -211,7 +211,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 x>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values
+            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
             p_eta_baro[j] = p_eta_baro[j]*np.exp(-x) + t * gj * c
 
 # Takes in a configuration of atoms (positions, velocities)

python_examples/md_nvt_lj.py

@@ -154,7 +154,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 x>0.001 else np.polyval([-1/24,1/6,-1/2,1.0],x) # Guard against small values
+            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
             p_eta[j] = p_eta[j]*np.exp(-x) + t * gj * c
 
 # Takes in a configuration of atoms (positions, velocities)