diff --git a/bd_nvt_lj.f90 b/bd_nvt_lj.f90 index 77d58d5..ca79d7c 100644 --- a/bd_nvt_lj.f90 +++ b/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 ) ) ) diff --git a/corfun.f90 b/corfun.f90 index 8f84b30..cf93c19 100644 --- a/corfun.f90 +++ b/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 diff --git a/diffusion_test.f90 b/diffusion_test.f90 index 2faf4cb..28def4b 100644 --- a/diffusion_test.f90 +++ b/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 diff --git a/error_calc.f90 b/error_calc.f90 index 4782c07..84504d1 100644 --- a/error_calc.f90 +++ b/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 diff --git a/md_npt_lj.f90 b/md_npt_lj.f90 index 7b0df55..a691a22 100644 --- a/md_npt_lj.f90 +++ b/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 diff --git a/md_nvt_lj.f90 b/md_nvt_lj.f90 index 6e51ec0..f0cc77c 100644 --- a/md_nvt_lj.f90 +++ b/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 diff --git a/python_examples/bd_nvt_lj.py b/python_examples/bd_nvt_lj.py index cac7db4..56bc62a 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 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) diff --git a/python_examples/corfun.py b/python_examples/corfun.py index 9b39a22..7317762 100755 --- a/python_examples/corfun.py +++ b/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) diff --git a/python_examples/diffusion_test.py b/python_examples/diffusion_test.py index a9ffde0..2f11fff 100755 --- a/python_examples/diffusion_test.py +++ b/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) diff --git a/python_examples/error_calc.py b/python_examples/error_calc.py index 266720f..8814f8b 100755 --- a/python_examples/error_calc.py +++ b/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) diff --git a/python_examples/md_npt_lj.py b/python_examples/md_npt_lj.py index 7667025..fb9c716 100755 --- a/python_examples/md_npt_lj.py +++ b/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) diff --git a/python_examples/md_nvt_lj.py b/python_examples/md_nvt_lj.py index 6d9577c..165f80b 100755 --- a/python_examples/md_nvt_lj.py +++ b/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)