From 1179380a38ce4b246401b622aeb9908352679ba9 Mon Sep 17 00:00:00 2001 From: Bruno Beltran Date: Thu, 18 Mar 2021 11:51:16 -0700 Subject: [PATCH] DOCS: more of ariana's suggestions --- doc/features.rst | 2 +- wlcsim/bd/__init__.py | 2 +- wlcsim/bd/homolog.py | 4 +- wlcsim/bd/rouse.py | 2 - wlcsim/bd/runge_kutta.py | 104 +++++++++++++++++++-------------------- 5 files changed, 53 insertions(+), 61 deletions(-) diff --git a/doc/features.rst b/doc/features.rst index 07f928d7..1caaaf73 100644 --- a/doc/features.rst +++ b/doc/features.rst @@ -3,7 +3,7 @@ Fortran Simulation Features ########################### -.. Why use wlcsim? Does it do what I need? +Why use wlcsim? Does it do what I need? Semiflexable polymers ===================== diff --git a/wlcsim/bd/__init__.py b/wlcsim/bd/__init__.py index 156052f4..bc23d8d7 100644 --- a/wlcsim/bd/__init__.py +++ b/wlcsim/bd/__init__.py @@ -4,7 +4,7 @@ A Brownian dynamics simulation requires specification of three things: 1. **The system state**. For us, this is :math:`\vec{r}_j(t_i) \in - \mathbb{r}^{N\times 3}`, the positions of the :math:`N` discrete beads + \mathbb{R}^{N\times 3}`, the positions of the :math:`N` discrete beads representing the polymer at each time :math:`t_i`. 2. **The forces**. Functions taking as input the system state and returning as output the forces on each bead. The main forces in any polymer diff --git a/wlcsim/bd/homolog.py b/wlcsim/bd/homolog.py index 00f5e46c..368976c3 100644 --- a/wlcsim/bd/homolog.py +++ b/wlcsim/bd/homolog.py @@ -3,8 +3,6 @@ See `homolog_points_to_loops_list` for a description of how the homologous network structure is handled. - -A specialized """ from ..plot import PolymerViewer from .rouse import recommended_dt @@ -156,7 +154,7 @@ def split_homologs_X(X, N, loop_list): return X1, X2 def split_homolog_x(x0, N, loop_list): - """Make an (N_tot,3) array from rouse_homologs into (2,N,3).""" + """Make an (N_tot,3) array from `rouse_homologs` into (2,N,3).""" N_tot, _ = x0.shape x1 = x0[:N,:].copy() x2 = np.zeros((N,3)) diff --git a/wlcsim/bd/rouse.py b/wlcsim/bd/rouse.py index 071781bc..0a7df9be 100644 --- a/wlcsim/bd/rouse.py +++ b/wlcsim/bd/rouse.py @@ -92,13 +92,11 @@ can be thought of as one large "effective" bead with a diffusivity :math:`N` times smaller than the individual beads, and will again exhibit regular diffusion (:math:`t^1` scaling). - It is worthwhile to note that the case of the continuous Rouse polymer, there is no shortest time scale (because of its fractal nature, wherein at each time point, the Rouse chain's spatial configuration traces out a Brownian walk). In this case, the :math:`t^{1/2}` behavior simply stretches out to arbitrarily short time scales. - For a more detailed discussion of the scaling behavior of Rouse MSD's, see the `measured_D_to_rouse` function. diff --git a/wlcsim/bd/runge_kutta.py b/wlcsim/bd/runge_kutta.py index 2e07420b..82038136 100644 --- a/wlcsim/bd/runge_kutta.py +++ b/wlcsim/bd/runge_kutta.py @@ -20,11 +20,20 @@ def rk4_thermal_lena(f, D, t, x0): - """x'(t) = f(x(t), t) + Xi(t), where Xi is thermal, diffusivity D. + r""" + Integrate using the method of Steph Weber (Phys. Rev. E., 82, 011913). + + Solves the equation + + .. math:: + + x'(t) = f(x(t), t) + \xi(t), - x0 is x(t[0]). + where :math:`\xi` represents thermal forces, and each dimension of ``x0`` + has diffusivity :math:`D`. - :math:`f: R^n x R -> R^n` + ``x0`` is :math:`x(t=0)`, and :math:`f: \mathbb{R}^n \times \mathbb{R} -> + \mathbb{R}^n` """ t = np.array(t) x0 = np.array(x0) @@ -53,51 +62,22 @@ def rk4_thermal_lena(f, D, t, x0): return x -def rk4_thermal_bruno(f, D, t, x0): +def euler_maruyama(f, D, t, x0): r""" - Test new method: Attempt to keep :math:`\omega` constant. + The most well-known stochastic integrator. - WARNING: does not converge strongly (autocorrelation function seems - higher than should be for OU process...), as is...x'(t) = f(x(t), t) + - Xi(t), where Xi is thermal, diffusivity D + Solves the equation - x0 is x(t[0]). + .. math:: - :math:`f: R^n x R -> R^n` - """ - t = np.array(t) - x0 = np.array(x0) - xsize = x0.shape - x = np.zeros(t.shape + x0.shape) - dts = np.diff(t) - x[0] = x0 - dxdt = np.zeros((4,) + x0.shape) # one for each RK step - Fbrown = np.sqrt(2*D / ((t[1] - t[0])/2)) - # at each step i, we use data (x,t)[i-1] to create (x,t)[i] - # in order to make it easy to pull into a new functin later, we'll call - # t[i-1] "t0", old x (x[i-1]) "x0", and t[i]-t[i-1] "h". - for i in range(1, len(t)): - h = dts[i-1] - t0 = t[i-1] - x0 = x[i-1] - x_est = x0 - dxdt[0] = f(x0, t0) + Fbrown # slope at beginning of time step - # random force estimate at midpoint - Fbrown = np.sqrt(2*D / ((t[i]-t[i-1])/2))*np.random.normal(size=xsize) - x_est = x0 + dxdt[0]*(h/2) # first estimate at midpoint - dxdt[1] = f(x_est, t0 + (h/2)) + Fbrown # estimated slope at midpoint - x_est = x0 + dxdt[1]*(h/2) # second estimate at midpoint - dxdt[2] = f(x_est, t0 + (h/2)) + Fbrown # second slope at midpoint - x_est = x0 + dxdt[2]*h # first estimate at next time point - # random force estimate at endpoint (and next start point) - Fbrown = np.sqrt(2*D / ((t[i]-t[i-1])/2))*np.random.normal(size=xsize) - dxdt[3] = f(x_est, t0 + h) + Fbrown # slope at end of time step - # final estimate is weighted average of slope estimates - x[i] = x0 + h*(dxdt[0] + 2*dxdt[1] + 2*dxdt[2] + dxdt[3])/6 - return x + x'(t) = f(x(t), t) + \xi(t), + where :math:`\xi` represents thermal forces, and each dimension of ``x0`` + has diffusivity :math:`D`. -def euler_maruyama(f, D, t, x0): + ``x0`` is :math:`x(t=0)`, and :math:`f: \mathbb{R}^n \times \mathbb{R} -> + \mathbb{R}^n` + """ t = np.array(t) x0 = np.array(x0) x = np.zeros(t.shape + x0.shape) @@ -116,7 +96,8 @@ def euler_maruyama(f, D, t, x0): def srk1_roberts(f, D, t, x0): - r"""From wiki, from A. J. Roberts. Modify the improved Euler scheme to + r""" + From wiki, from A. J. Roberts. Modify the improved Euler scheme to integrate stochastic differential equations. [1], Oct 2012. If we have an Ito SDE given by @@ -129,13 +110,22 @@ def srk1_roberts(f, D, t, x0): .. math:: - \vec{K}_1 = h \vec{a}(t_k, \vec{X}_k) + (\Delta W_k - S_k\sqrt{h}) \vec{b}(t_k, \vec{X}_k) - \vec{K}_2 = h \vec{a}(t_{k+1}, \vec{X}_k + \vec{K}_1) + (\Delta W_k - S_k\sqrt{h}) \vec{b}(t_{k+1}, \vec{X}_k + \vec{K}_1) + \vec{K}_1 = h \vec{a}(t_k, \vec{X}_k) + (\Delta W_k - S_k\sqrt{h}) + \vec{b}(t_k, \vec{X}_k) + + .. math:: + + \vec{K}_2 = h \vec{a}(t_{k+1}, \vec{X}_k + \vec{K}_1) + + (\Delta W_k - S_k\sqrt{h}) \vec{b}(t_{k+1},\vec{X}_k + \vec{K}_1) + + .. math:: + \vec{X}_{k+1} = \vec{X}_k + \frac{1}{2}(\vec{K}_1 + \vec{K}_2) where :math:`\Delta W_k = \sqrt{h} Z_k` for a normal random :math:`Z_k \sim N(0,1)`, and :math:`S_k=\pm1`, with the sign chosen uniformly at random - each time.""" + each time. + """ t = np.array(t) x0 = np.array(x0) x = np.zeros(t.shape + x0.shape) @@ -164,8 +154,14 @@ def srk1_roberts(f, D, t, x0): # simple test case def ou(x0, t, k_over_xi, D, method=rk4_thermal_lena): - "simulate ornstein uhlenbeck process with theta=k_over_xi and sigma^2/2=D" - def f(x,t): + r""" + Test simulation: the Ornstein-Uhlenbeck process. + + Simulates a standard Ornstein-Uhlenbeck process with ``theta=k_over_xi`` + and ``sigma^2/2=D``. This is the appropriate parameterization to match the + limiting case of ``N=2`` in `wlcsim.bd.rouse.with_integrator`. + """ + def f(x, t): return -k_over_xi*x return method(f, D=D, t=t, x0=x0) @@ -179,7 +175,7 @@ def _get_scalar_corr(X): for i in range(num_samples): for j in range(num_t): for k in range(j, num_t): - corr[k-j] = corr[k-j] + X[i,k]*X[i,j] + corr[k-j] = corr[k-j] + X[i, k]*X[i, j] count[k-j] = count[k-j] + 1 return corr, count @@ -193,7 +189,7 @@ def _get_vector_corr(X): for i in range(num_samples): for j in range(num_t): for k in range(j, num_t): - corr[k-j] = corr[k-j] + X[i,k]@X[i,j] + corr[k-j] = corr[k-j] + X[i, k]@X[i, j] count[k-j] = count[k-j] + 1 return corr, count @@ -211,7 +207,7 @@ def _get_bead_msd(X, k=None): count = np.zeros((num_t,)) for i in range(num_t): for j in range(i, num_t): - ta_msd[j-i] += (X[j,k] - X[i,k])@(X[j,k] - X[i,k]) + ta_msd[j-i] += (X[j, k] - X[i, k])@(X[j, k] - X[i, k]) count[j-i] += 1 return ta_msd, count @@ -219,15 +215,15 @@ def _get_bead_msd(X, k=None): @jit(nopython=True) def _msd(x): result = np.zeros_like(x) - for delta in range(1,len(x)): - for i in range(delta,len(x)): + for delta in range(1, len(x)): + for i in range(delta, len(x)): result[delta] += (x[i] - x[i-delta])**2 result[delta] = result[delta] / (len(x) - delta) return result -# test different integrators below on simply OU process def test_ou_autocorr(method=srk1_roberts): + """Test an individual integrator visually using an OU process.""" import scipy.stats import matplotlib.pyplot as plt k = 2