wlcsim Python Module

Compute results related to (twistable, stretchable, shearable) wormlike chains.

This package can be largely separated into

1. Modules for processing simulation output from our FORTRAN codebase (wlcsim.data, wlcsim.input) 2. Modules for generating polymer chains at equilibrium (both directly in wlcsim.chains and using Monte Carlo wlcsim.mc). 3. Modules for simulating diffusing polymers (wlcsim.bd). 4. Modules for plotting these various chains and their properties (wlcsim.plot and wlcsim.plot_wlcsim).

wlcsim.bd

wlcsim.bd.rouse

Simulate Rouse polymers.

Notes

There are various ways to parameterize Rouse polymers. In this module, we use the convention that N is the number of beads of a Rouse polymer (contrary to e.g. Doi & Edwards, where N is the number of Kuhn lengths in the polymer). We instead use Nhat for the number of Kuhn lengths in the polymer.

D is the “diffusivity of a Kuhn length”, i.e. kbT/xi, where xi is the dynamic viscosity of the medium in question, as typically found in the Langevin equation for a Rouse polymer

\[\xi \frac{d}{dt} \vec{r}(n, t) = k \frac{d^2}{dn^2} \vec{r}(n, t) + f^{(B)}\]

This means that to simulate a Rouse polymer diffusing in a medium with viscosity xi, the diffusion coefficient of each bead should be set to Dhat = D (N/Nhat). Since \(\xi\) is in units of “viscosity per Kuhn length” and (N/Nhat) is in units of “number of beads over number of Kuhn lengths”, this can be thought of as changing units from “viscosity per Kuhn length” to “viscosity per bead”.

Some books use b for the mean (squared) bond length between beads in the discrete gaussian chain. We instead use b to be the real Kuhn length of the polymer, so that the mean squared bond length bhat**2 is instead given by L0*b, where L0 is L0 = L/(N-1) is the amount of polymer “length” represented by the space between two beads.

The units chosen here make the most sense if you don’t actually consider the polymer to be a “real” Rouse polymer, but instead an e.g. semiflexible chain with a real “length” whose long distance statistics are being captured by the simulation.

To compare these results to the rouse.analytical.rouse module, use

>>> plt.plot(t_save, sim_msd)
>>> plt.plot(t_save, wlcsim.analytical.rouse.rouse_mid_msd(t_save, b, Nhat, D,
...                                                        num_modes=int(N/2)))
>>> plt.plot(t_save, 6*Dhat*t_save)
>>> plt.plot(t_save, 6*(Dhat/N)*t_save)  # or 6*(D/Nhat)*t_save
>>> # constant prefactor determined empirically...
>>> # no, using the value of 6*(Dhat/N)*t_R doesn't work...
>>> plt.plot(t_save, 1.9544100*b*np.sqrt(D)*np.sqrt(t_save))

where cutting off the number of modes corresponds to ensuring that the rouse behavior only continues down to the length scale of a single “bead”, thus matching the simulation at arbitrarily short time scales. (Notice that we abide by the warning above. Our Nhat is wlcsim.analytical.rouse’s N).

Example

For example, if you want to simulate a megabase of DNA, which has a real linear length of about 1e6 megabase/base * 0.33 nm/base ~ 3e5 nm, and a Kuhn length of 1e2 nm, then

wlcsim.bd.rouse.jit_confined_srk1

Add an elliptical confinement.

Energy is like cubed of distance outside of ellipsoid, pointing normally back in. times some factor Aex controlling its strength.

wlcsim.bd.rouse.jit_confinement_clean

Unfortunately, by “cleaning up” this function, we also make it 2x slower.

#worthit

wlcsim.bd.rouse.jit_srk1

Faster version of wlcsim.bd.rouse.rouse using jit.

N=101,L=100,b=1,D=1 takes about 3.5min to run when t=np.linspace(0, 1e5, 1e7+1) adding t_save does not slow function down

Our srk1 scheme is accurate as long as \(\Delta t\) is less than the transition time where the MSD goes from high k behavior (\(t^1\)) to Rouse behavior (\(t^{1/2}\)). This is exactly the time required to diffuse a Kuhn length, so we just need \(\Delta t < \frac{\hat{b}^2}{6\hat{D}}\) in 3D. the “crossover” from fast-k to rouse-like behavior takes about one order of magnitude in time, but adding more than that doesn’t seem to make the MSD any more accurate, so we suggest setting \(\Delta t\) to one tenth of the bound above.

the number of orders of magnitude of “rouse” scaling the simulation will capture is exactly dictated by the ratio between this time scale and the rouse relaxation time (so like N**2?)

recall (doi & edwards, eq 4.25) that the first mode’s relaxation time is \(\tau_1 = \frac{\xi N^2}{k \pi^2 }\). and the \(p\)th mode is \(\tau_p = \tau_1/p^2\) (this is the exponential falloff rate of the \(p\)th mode’s correlation function).

wlcsim.bd.rouse.measured_D_to_rouse(Dapp, d, N, bhat=None, regime='rouse')

Get the full-polymer diffusion coefficient from the “apparent” D.

In general, a discrete Rouse polymer’s MSD will have three regimes. On time scales long enough that the whole polymer diffuses as a large effective particle, the MSD is of the form \(6 D/\hat{N} t\), where, as is true throughout this module, \(D\) is the diffusivity of a Kuhn length, and \(\hat{N}\) is the number of Kuhn lengths in the polymer.

This is true down to the full chain’s relaxation time (the relaxation time of the 0th Rouse mode, also known as the “Rouse time”) \(t_R = N^2 b^2 / D\). For times shorter than this, the MSD will scale as \(t^{1/2}\). Imposing continuity of the MSD, this means that the MSD will behave as \(\kappa_0 D/\hat{N} (t_R t)^{1/2}\), where \(\kappa_0\) is a constant that I’m too lazy to compute using \(\lim_{t\to 0}\) of the analytical Rouse MSD result, so I just determine it empirically. We rewrite this MSD as \(\kappa b D^{-1/2} t^{1/2}\), and find by comparing to the analytical Rouse theory that \(\kappa = 1.9544100(4)\).

Eventually (at extremely short times), most “real” polymers will eventually revert to a MSD that scales as \(t^1\). This cross-over time/length scale defines a “number of segments” \(\tilde{N}\), where the diffusivity of a segment of length \(\tilde{L} = L/(\tilde{N} - 1)\) matches the time it takes stress to propagate a distance \(\tilde{L}\) along the polymer. Said in other words, for polymer lengths smaller than \(\tilde{L}\), the diffusivity of the segment outruns the stress communication time between two neighboring segments.

The diffusivity at these extremely short times will look like \(6 D (\tilde{N} / \hat{N}) t\). In order to simulate this behavior exactly, one can simply use a polymer with \(\tilde{N}\) beads, then all three length scales of the MSD behavior will match.

This three-regime MSD behavior can be very easily visualized in log-log space, where it is simply the continuous function defined by the following three lines. From shortest to longest time scales, \(\log{t} + \log{6D(\tilde{N}/\hat{N})}\), \((1/2)\log{t} + \log{\kappa b D^{1/2}}\), and \(\log{t} + \log{6D\hat{N}}\). The lines (in log-log space), have slopes 1, 1/2, and 1, respectively, and their “offsets” (y-intercept terms with the log removed) are typically referred to as \(D_\text{app}\).

The simulations in this module use the diffusivity of a single Kuhn length as input (i.e. plain old \(D\) is expected as input) but typically, measurements of diffusing polymers are done below the Rouse time, and for long enough polymers (like a mammalian chromosome), it may be difficult to impossible to measure unconfined diffusion at time scales beyond the Rouse time. This means you often can’t just look at the diffusivity of the whole polymer and multiply by \(\hat{N}\) to get the diffusivity of a Kuhn length. And since there’s not really a principled way in general to guess \(\tilde{N}\), you usually can’t use the short-time behavior to get \(D\) either, even supposing you manage to measure short enough times to see the cross-over back to \(t^1\) scaling. This means that usually the only way one can extract \(D\) is by measuring \(D_\text{app}\), since we can use the fact that \(D_\text{app} = \kappa b D^{1/2}\) (the \(\kappa\) value quoted above only works for 3-d motion. I haven’t bothered to check if it scales linearly with the number of dimensions or like \(\sqrt{d}\) for \(d\)-dimensional motion, but this shouldn’t be too hard if you are interested).

In short, this function gives you \(D\) given \(D_\text{app}\).

Parameters:

• D_app (float) – The measured \(D_\text{app} based on the y-intercept of the MSD in log-log space. This can be computed as \), assuming that t is purely in the regime where the MSD scales like \(t^{1/2}\).
• Nhat (int) – number of beads in the polymer
• d (int) – number of dimensions of the diffusion (e.g. 3 in 3-d space)

Returns:

D_rouse – the diffusivity of a Kuhn length of the polymer

Return type:

float

Notes

What follows is an alternate way of “thinking” about the three regimes of the MSD, from a simulation-centric perspective, that I typed up back when I first wrote this function, but has been obsoleted by the description at the start of this docstring.

In general, a discrete Rouse polymer’s MSD will have three regimes. On time scales short enough that chain connectivity is not relevant, the MSD is of the form \(6 \hat{D} t\).

As soon as chain connectivity begins to affect the dynamics, the form of the MSD will change to \(\kappa \hat{b} \hat{D}^{1/2} t^{1/2}\), where we have determined the constant \(\kappa\) empirically (using our exact analytical theory) to be approximately 1.9544100(4). Note that (up to a constant factor), this is just \(6 \hat{D} (t_R t)^{1/2}\), where \(t_R\) is the relaxation time (aka Rouse time) of the polymer, given by \(\hat{b}^2 \tilde{N}^2 / \hat{D}\).

This is because the Rouse behavior only continues up to the relaxation time \(t = t_R\), where the MSD transitions into the form \(6 \frac{\hat{D}}{\tilde{N}} t\). If you ask what “line” with slope 1/2 in log-log space intersects \(6 \frac{\hat{D}}{\tilde{N}} t\) at \(t = t_R\), you get the above form for the sub-Rouse time MSD.

Here, \(\tilde{N}\) is distinguished from \(N\) only as a formality. For a simulation of a discrete Rouse polymer with a fixed number of beads, \(N = \tilde{N}\) is exactly the number of beads. For a “real” polymer, it is simply a numerical parameter that describes how long the Rouse behavior lasts (in the limiit of short times). For a “true” Rouse polymer (the fractal object) N is infinity and there is no time scale on which the MSD transitions back into \(t^1\) behavior. But, for any real polymer, this will “almost always” eventually happen.

wlcsim.bd.rouse.recommended_dt(N, L, b, D)

Recommended “dt” for use with rouse*jit family of functions.

Currently set to \(\frac{1}{10}\frac{b^2}{6D}\).

See the jit_srk1() docstring for source of this time scale.

wlcsim.bd.rouse.with_integrator(N, L, b, D, t, x0=None, integrator=<function rk4_thermal_lena>)

Simulate a Rouse polymer made of N beads free in solution.

Parameters:

• N (float) – Number of beads in the chain.
• L (float) – Length of chain.
• b (float) – Kuhn length of the chain (same units as L).
• D (float) – Diffusion coefficient. (Units of length**2/time). In order to compute D from a single-locus MSD, use measured_D_to_rouse.
• t ((Nt,) array_like of float) – Time points to use for stepping the integrator. None of our current integrators is currently capable of taking time steps larger than the time scale of the highest Rouse mode of the finite chain (use ~.bd.recommended_dt to compute the optimal dt).
• x0 ((N, 3) array_like of float, optional) – The initial locations of each bead. If not specified, defaults to initializing from the free-draining equilibrium, with the first bead at the origin.
• integrator (Callable[[ForceFunc, float, Times, Positions], Positions]) – Either ~.runge_kutta.rk4_thermal_lena or ~.runge_kutta.srk1_roberts

Returns:

The positions of the N beads at each of the Nt time points.

Return type:

(Nt, N, 3) array_like of float

Notes

The polymer beads can be described by the discrete Rouse equations

\[\xi \frac{dx(i, t)}{dt} = - k (x(i, t) - x(i+1, t)) - k (x(i, t) - x(i-1, t)) + R(t)\]

where \(\xi = k_BT/D\), \(k = 3k_BT/b^2\), \(b\) is the Kuhn length of the polymer, \(D\) is the self-diffusion coefficient of a bead, and \(R(t)/\xi\) is a delta-correlated stationary Gaussian process with mean zero and \(\langle R(t) R(t')/\xi^2 \rangle = 2DI\delta(t-t')\).

Notice that in practice, \(k/\xi = 3D/b^2\), so we do not need to include mass units (i.e. there’s no dependence on \(k_BT\)).

wlcsim.bd.homolog

BD of two homologously-linked linear Rouse polymers.

See homolog_points_to_loops_list for a description of how the homologous network structure is handled.

A specialized

wlcsim.bd.homolog._jit_rouse_homologs

“Inner loop” of rouse_homologs.

Some typing stuff: needs N, N_tot as int, needs tether_list/loop_list as dtype == int64, and uses loop_list.shape[0] to get num_loops. This is in particular important for empty list cases, where you need to do stuff like

>>> x = rouse._jit_rouse_homologs(int(N), int(2*N),
>>>         np.array([]).astype(int),
>>>         np.array([[]]).T.astype(int),
>>>         L, b, D, Aex, R, R, R, t, t_save)

wlcsim.bd.homolog.points_to_loops_list(N, homolog_points)

Create loops list for jit_rouse_linked.

Parameters:

• N (int) – number of beads in each polymer if they were unbound
• homolog_points ((L,) array_like) – list of loci indices \({l_i}_1^L, l_i\in[1,N]\) (beads) that are homologously paired.

Returns:

• N_tot (int) – final size of the x0 array for the simulation.
• loop_list ((L,4) np.array) – data structure used by jit_rouse_linked to correctly compute the spring forces in the network.

Notes

You can use the output like

>>> k1l, k1r, k2l, k2r = loop_list[i]

where “k[i][l,r]” refers to the bead demarcating the [l]eft or [r]ight end of a given “homologous loop” on each of the polymers polymer \(i\in[1,2]\).

Namely, the connectivity structure of our polymer “network” is that for each k1l, k1r, k2l, k2r in loop_list, the beads k1l thru k1r are linked sequentially (as are k2l thru k2r), but the beads k1[l,r] are linked to the beads k2[l,r], respectively. This network structure can represent arbitrary “homologously linked” polymers, i.e. those of the form

--  ---  -------  --  -
  \/   \/       \/  \/
  /\   /\       /\  /\
--  ---  -------  --  -

Notice in the example below that loop_list[1:,1] and loop_list[:-1,0] are both the same as homolog_list.

Example

Here’s an example for a typical case. With 101 beads per polymer (so 100 inter-bead segments), and one tenth of the beads bound to each other, we might get the following Poisson-space homolog points:

>>> N = 101; FP = 0.1

>>> homolog_points = np.where(np.random.rand(int(N)) < FP)[0]

>>> homolog_points
array([ 5, 18, 27, 59, 68, 82, 90, 94])

These eight homologous junctions would lead to a simulation array length of 2*N - len(homolog_points) = 202 - 8 = 194. So N_tot will be 194, and the loop_list will look like

>>> N_tot, loop_list = rouse.homolog_points_to_loops_list(N, homolog_points)

>>> loop_list
array([[ -1,   5, 101, 105],
        [  5,  18, 106, 117],
        [ 18,  27, 118, 125],
        [ 27,  59, 126, 156],
        [ 59,  68, 157, 164],
        [ 68,  82, 165, 177],
        [ 82,  90, 178, 184],
        [ 90,  94, 185, 187],
        [ 94, 101, 188, 193]])

This can be read as follows. The first four beads of each polymer are “free”, so first row is saying that beads 0 thru 4 correspond to unlinked beads in polymer 1. bead 5 is the linked bead (has half the diffusivity). since the polymer is of length 101, if there were no connections, the first polymer would extend from bead 0 to bead 100, and the second from bead 101 to bead 201. but since bead “5” is linked to (what would be) bead 106, that bead is omitted from the second polymer, and beads 101-105 are the “free end” of the second polymer.

The second row corresponds to the first actual “loop”. The two tethering points are beads 5 and 18 of the first polymer. What would have been beads 106 and 119 (that’s beads 5 and 18 of the second polymer) are omitted from the list of simulated beads. So bead 106 is linked to bead 5 and bead 117 is linked to bead 18, and beads 106-117 are linked sequentially.

Similarly, beads 18 and 118 (and 27 and 125) are linked, with beads 118-125 being linked sequentially.

wlcsim.bd.homolog.rouse(N, FP, L, b, D, Aex, rx, ry, rz, t, t_save=None, tether_list=None)

BD simulation of two homologous yeast chromosomes in meiosis.

An arbitrary elliptical-shaped confinement can be chosen, and in order to simulate synaptonemal formation in prophase, some fraction FP of the “homologous” beads of the polymer can be “rigidly tethered” to each other. The simulation treats these pairs as being one bead (with half the diffusion coefficient). The elastic force can then be communicated between the two polymers via the connecting loci.

Depending on what part of prophase is being simulated, either (or both) of the telomeres or centromeres can be tethered to the confinement by including them in tether_list. Currently the confinement force and the tethering force are both cubic w.r.t. the distance from the confining ellipse with the same strength (Aex), and have the form described in f_conf_ellipse() (and f_tether_ellipse()).

Parameters:

• N (int) – Number of beads to use in the discretized Rouse chain
• FP (float) – FP\(\in[0,1] is fraction of beads in the chain (out of \)) that will be tethered to each other
• L (float) – Length of the chain (in desired length units)
• b (float) – Kuhn length of the chain (in same length units as L)
• D (float) – The diffusivity of a Kuhn length of the chain. This can often difficult to measure directly, but see the documentation for the function measured_D_to_rouse() for instructions on how to compute this value given clean measurements of short-time regular diffusion or Rouse-like MSD behavior
• Aex (float) – multiplicative prefactor controlling strength of confinment and tethering forces
• ry, rz (rx,) – radius of confinement in x, y, z direction(s) respectively
• t ((M,) float, array_like) – The BD integrator will step time from t[0] to t[1] to t[2] to …
• t_save ((m,) float, array_like) – t_save\(\subset\) will be the time steps where the simulation output is saved. Default is t_save = t.
• tether_list (int, array_like) – list of bead indices that will be tethered to the confinement surface, for now, both polymers must be tethered at the same beads (but this would not be hard to change).

Returns:

• tether_list (List<int>) – Indices (into compressed Ntot output) that were tethered to the confinement surface
• loop_list ((~N*FP, 4) int) – Output of points_to_loops_list used to specify which points are tethered
• X ((len(t_save), Ntot, 3)) – output of BD simulation at each time in t_save. The number of output beads Ntot <= 2*N excludes the beads on the second polymer whose positions are determined by the fact that they are tethered to the first polymer. Use split_homologs_X() to reshape this output into shape (2, len(t_save), N, 3).

wlcsim.bd.homolog.sim_loc(i, N, loop_list)

Return actual indices for bead “i” of the second polymer.

Useful for indexing directly into the trucated simulation output array.

Parameters:

• i ((M,) int, array_like) – beads index is requested for
• N (int) – actual number of beads in each polymer
• loop_list ((~N*FP, 4) int) – connectivity of the homologous polymers being simulated

Returns:

iloc – index in simulation array (or None if the bead is a paired bead) for each bead in i.

Return type:

(M,) Optional<int>, array_like

wlcsim.bd.homolog.split_homolog_x(x0, N, loop_list)

Make an (N_tot,3) array from rouse_homologs into (2,N,3).

wlcsim.bd.homolog.split_homologs_X(X, N, loop_list)

Make the (Nt,Ntot,3) output of rouse_homologs to (2,Nt,N,3).

wlcsim.bd.forces

Force fields used in our BD integrators.

Our BD code can be summarized as typically looking like:

X = init_function()
for t_i in requested_times:
    X = integrator_step(f_all, X, prev_t, t_i)
    prev_t = t_i

where f_all takes X (size N_beads by N_dimensions) and t and returns the total force on each of the beads at time t.

This module holds all the functions that can be used to construct f_all. Broadly, they are named based on the source of the force:

• f_elas_* are elastic forces, they define the physics of the polymer being simulated. Linear, ring, and different types of network polymers have their own independent implementations.
• f_conf_* are confinement forces.
• f_tether_* tether individual beads of the polymer to an external structure.
• Others may be added in the future.

Notes

Because they are in the “innermost loop”, it is important that these functions are just-in-time compile-able when possible.

wlcsim.bd.forces.f_conf_ellipse

Compute soft (cubic) force due to elliptical confinement.

wlcsim.bd.forces.f_elas_homolog_rouse

Compute spring forces on two “homologously linked linear rouse polymers.

The two polymers are (homologously) at beads specified by loop_list. While each polymer is of length N beads, the beads that are hooked together move together identically, so you have x0.shape[0] == 2*N - len(loop_list)

We assume that, per the spec in ~.homolog_points_to_loops_list, a polymer of the shape

++   ^^^   '''''''   ""   =
  \ /   \ /       \ /  \ /
   .     .         .    .
  / \   / \       / \  / \
--   ---   -------   --   -

will be laid out in memory like

--.---.-------.--.-++^^^'''''''""=

If we call the lower polymer “polymer 1”, and all the beads on the top row “polymer 2”, then it makes sense to first compute the spring forces between adjacent beads in polymer 1 as if there was no “polymer 2”. Then, we can compute all the forces between beads in “polymer 2” and the “.” beads. Finally, for each “loop” (represented as symbols of different types above) in polymer 2, we can compute the usual linear forces along the polymer.

wlcsim.bd.forces.f_elas_linear_rouse

Compute spring forces on single, linear rouse polymer.

wlcsim.bd.forces.f_elas_linear_sswlc

Compute spring forces and torques on each bead of dsswlc.

wlcsim.bd.forces.f_tether_ellipse

Compute soft (cubic) force towards an elliptical confinement.

wlcsim.bd.runge_kutta

Some spare stochastic integrators.

Not currently used because if you pass a function as an argument you can’t @numba.jit that.

For rouse chains, Bruno (that’s me) has tested the srk1 integrator and found it to work extremely well (see notes in wlcsim.bd.rouse on suggested dt).

For WLC/ssWLC, Lena found it useful to use a scheme that was higher deterministic order in time, to resolve the high elastic energies involved. Bruno has not yet tested this explicitly.

You can find Lena’s algorithm (hold brownian force constant over an rk4 step) below. It is unlikely to be strongly convergent except if you subsample below the actual desired time resolution.

wlcsim.bd.runge_kutta._get_bead_msd

center bead by default

for 1e4-long time arrays, this takes ~10-30s on my laptop

wlcsim.bd.runge_kutta._get_scalar_corr

fast correlation calculation for testing

wlcsim.bd.runge_kutta._get_vector_corr

fast correlation calculation for testing

wlcsim.bd.runge_kutta.ou(x0, t, k_over_xi, D, method=<function rk4_thermal_lena>)

simulate ornstein uhlenbeck process with theta=k_over_xi and sigma^2/2=D

wlcsim.bd.runge_kutta.rk4_thermal_bruno(f, D, t, x0)

Test new method: Attempt to keep \(\omega\) constant.

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

x0 is x(t[0]).

\(f: R^n x R -> R^n\)

wlcsim.bd.runge_kutta.rk4_thermal_lena(f, D, t, x0)

x’(t) = f(x(t), t) + Xi(t), where Xi is thermal, diffusivity D.

x0 is x(t[0]).

\(f: R^n x R -> R^n\)

wlcsim.bd.runge_kutta.srk1_roberts(f, D, t, x0)

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

\[d\vec{X} = \vec{a}(t, \vec{X}) + \vec{b}(t, \vec{X}) dW\]

then

\[\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{X}_{k+1} = \vec{X}_k + \frac{1}{2}(\vec{K}_1 + \vec{K}_2)\]

where \(\Delta W_k = \sqrt{h} Z_k\) for a normal random \(Z_k \sim N(0,1)\), and \(S_k=\pm1\), with the sign chosen uniformly at random each time.

wlcsim.bd.tsswlc

Simulation twistable, stretchable-shearable wormlike chains.

Notes

The ssWLC is implemented as described in Koslover et. al., Soft Matter, 2013, 9, 7016. This defines the locations (\(r^{(i)}(t_j)\)) and tangent vectors (\(u^{(i)}(t_j)\)) of each (\(i\)th) bead at each time point (\(t_j\)).

The twist energy is (loosely speaking) quadratic in the angle that one must rotate the material normals of one bead to get them to align with the next bead (after the two beads have been rotated once to make the tangent vectors align).

There are other strategies one might take for this, but we will leave them as “#TODO” in the code for now.

wlcsim.bd.tsswlc._jit_sswlc_clean(N, L, lp, t, t_save, e_b, gam, e_par, e_perp, eta, xi_u, xi_r)

WARNING: known to be buggy…

wlcsim.bd.tsswlc._jit_sswlc_lena(N, L, lp, t, t_save, e_b, gam, e_par, e_perp, eta, xi_u, xi_r)

TODO: test

wlcsim.bd.tsswlc.sswlc(N, L, lp, t, t_save)

Simulate a stretchable-shearable WLC.

Parameters:

• N (int) – The number of beads to use.
• L (float) – The length of the chain.
• lp (float) – The persistence length of the underlying WLC being approximated. (Must be same units as L.
• t (array_like) – The times for which to simulate the polymer’s motion.
• t_save ((M, ) array_like) – The subset of t for which we should save the output.

Returns:

• r ((M, N, 3) array of float) – The positions of each bead at each save time.
• u ((M, N, 3) array of float) – The tangent vectors to the underlying wormlike chain at each bead.

wlcsim.analytical

Modules containing routines to evaluate analytical results for various polymer chains (equilibrium statistics and dynamics).

Rouse polymer, analytical results.

Notes

There are two parameterizations of the “Rouse” polymer that are commonly used, and they use the same variable name for two different things.

In one, N is the number of Kuhn lengths, and in the other, N is the number of beads, each of which can represent an arbitrary number of Kuhn lengths.

wlcsim.analytical.rouse.end2end_distance(r, lp, N, L)

For now, always returns values for r = np.linspace(0, 1, 50001).

Parameters:

• r ((N, ) float, array_like) – the values at which to evaluate the end-to-end probability distribution. ignored for now (TODO: fix)
• lp (float) – persistence length of polymer
• N (int) – number of beads in polymer
• L (float) – polymer length

Returns:

• x ((5001,) float) – np.linspace(0, 1, 5001)
• g ((5001,) float) – \(P(|R| = x | lp, N, L)\)

Notes

Uses the gaussian chain whenever applicable, ssWLC tabulated values otherwise. If you request parameters that require a WLC end-to-end distance, the function will ValueError.

wlcsim.analytical.rouse.end2end_distance_gauss(r, b, N, L)

in each dimension… ? seems to be off by a factor of 3 from the simulation….

wlcsim.analytical.rouse.end_to_end_corr(t, D, N, num_modes=10000)

Doi and Edwards, Eq. 4.35

wlcsim.analytical.rouse.gaussian_G

Green’s function of a Gaussian chain at N Kuhn lengths of separation, given a Kuhn length of b

wlcsim.analytical.rouse.gaussian_Ploop

Looping probability for two loci on a Gaussian chain N kuhn lengths apart, when the Kuhn length is b, and the capture radius is a

wlcsim.analytical.rouse.kp_over_kbt

“non-dimensionalized” k_p is all that’s needed for most formulas, e.g. MSD.

Type:k_p/(k_B T)

wlcsim.analytical.rouse.linear_mid_msd

modified from Weber Phys Rev E 2010, Eq. 24.

wlcsim.analytical.rouse.linear_mscd(t, D, Ndel, N, b=1, num_modes=10000)

Compute mscd for two points on a linear polymer.

Parameters:

• t ((M,) float, array_like) – Times at which to evaluate the MSCD
• D (float) – Diffusion coefficient, (in desired output length units). Equal to \(k_BT/\xi\) for \(\xi\) in units of “per Kuhn length”.
• Ndel (float) – Distance from the last linkage site to the measured site. This ends up being (1/2)*separation between the loci (in Kuhn lengths).
• N (float) – The full lengh of the linear polymer (in Kuhn lengths).
• b (float) – The Kuhn length (in desired length units).
• num_modes (int) – how many Rouse modes to include in the sum

Returns:

mscd – result

Return type:

(M,) np.array<float>

wlcsim.analytical.rouse.ring_mscd(t, D, Ndel, N, b=1, num_modes=10000)

Compute mscd for two points on a ring.

Parameters:

• t ((M,) float, array_like) – Times at which to evaluate the MSCD.
• D (float) – Diffusion coefficient, (in desired output length units). Equal to \(k_BT/\xi\) for \(\xi\) in units of “per Kuhn length”.
• Ndel (float) – (1/2)*separation between the loci on loop (in Kuhn lengths)
• N (float) – full length of the loop (in Kuhn lengths)
• b (float) – The Kuhn length, in desired output length units.
• num_modes (int) – How many Rouse modes to include in the sum.

Returns:

mscd – result

Return type:

(M,) np.array<float>

wlcsim.analytical.rouse.rouse_large_cvv_g(t, delta, deltaN, b, D)

Cvv^delta(t) for infinite polymer.

Lampo, BPJ, 2016 Eq. 16.

wlcsim.analytical.rouse.rouse_mode

Eigenbasis for Rouse model.

Indexed by p, depends only on position n/N along the polymer of length N. N=1 by default.

Weber, Phys Rev E, 2010 (Eq 14)

wlcsim.analytical.rouse.rouse_mode_coef

Weber Phys Rev E 2010, after Eq. 18.

Type:k_p

wlcsim.analytical.rouse.test_rouse_msd_line_approx()

Figure out what Dapp is exactly using our analytical result.

if we use msd_approx(t) = 3*bhat*np.sqrt(Dhat*t)/np.sqrt(3)*1.1283791(6) then |msd(t) - msd_approx(t)|/msd(t) = np.sqrt(2)/N*np.power(t, -1/2). msd_approx(t) > msd(t) in this range, so that means that msd_approx(t)/(np.sqrt(2)/N*np.power(t, -1/2) + 1) = msd(t).

if we redefine msd_approx with this correction, the new relative error is about np.sqrt(2)/N/100*t**(-1/2)? gotta see if this carries over to other chain parameters though…it does not… this time msd(t) is bigger, so msd_approx = msd_approx(t)/(1 - np.sqrt(2)/N/100*t**(-1/2)).

okay actually looks like extra factor is additive?

for N = 1e8+1; L = 25; b = 2; D = 1;, we have msd_approx(t) - msd(t) = 1.01321183e-07

for N = 1e8+1; L = 174; b = 150; D = 166;, we have msd_approx(t) - msd(t) = 5.2889657(3)e-05

for N = 1e8+1; L = 174; b = 15; D = 166;, we have msd_approx(t) - msd(t) = 5.2889657(3)e-06

for N = 1e8+1; L = 17.4; b = 15; D = 166;, we have msd_approx(t) - msd(t) = 5.2889657(3)e-07

for N = 1e8+1; L = 17.4; b = 15; D = 16.6;, we have msd_approx(t) - msd(t) = 5.2889657(3)e-07 (no change)

for N = 1e7+1; L = 17.4; b = 15; D = 16.6;, we have msd_approx(t) - msd(t) = 5.2889657(3)e-06 (no change)

so the answer is like 3*bhat*np.sqrt(Dhat*t)/np.sqrt(3)*1.1283791615 - 0.202642385398*b*L/N

Fractional (Rouse) polymer. A Rouse polymer in a viscoelastic medium.

For the case alpha=1, these formulas should reduce to the analogous ones in wlcsim.analytical.rouse.

wlcsim.analytical.fractional.frac_cv(t, alpha, D)

Velocity autocorrelation of a fractionally-diffusing particle. Weber, Phys Rev E, 2010 (Eq 32)

wlcsim.analytical.fractional.frac_discrete_cv(t, delta, alpha, D)

Discrete velocity autocorrelation of a fractionally-diffusing particle. Weber, Phys Rev E, 2010 (Eq 33)

wlcsim.analytical.fractional.frac_discrete_cv_normalized(t, delta, alpha)

Normalized discrete velocity autocorrelation of a fractionally-diffusing particle. Should be equivalent to

frac_discrete_cv(t, delta, 1, 1)/frac_discrete_cv(0, delta, 1, 1)

Lampo, BPJ, 2016 (Eq 5)

wlcsim.analytical.fractional.frac_msd(t, alpha, D)

MSD of fractionally diffusing free particle.

Weber, Phys Rev E, 2010 (Eq 10)

wlcsim.analytical.fractional.frac_rouse_cvv(t, delta, n1, n2, alpha, b, N, D, min_modes=500, rtol=1e-05, atol=1e-08, force_convergence=True)

Velocity cross-correlation of two points on fractional Rouse polymer

rtol/atol specify when to stop adding rouse modes. a particle (t,delta) pair is considered to have converged when np.isclose returns true give rtol/atol and p is even (p odd contributes almost nothing)

Lampo, BPJ, 2016 Eq. 10.

wlcsim.analytical.fractional.frac_rouse_mid_msd(t, alpha, b, N, D, num_modes=1000)

Weber Phys Rev E 2010, Eq. 24.

wlcsim.analytical.fractional.frac_rouse_mode_corr(p, t, alpha, b, N, D)

Weber Phys Rev E 2010, Eq. 21.

wlcsim.analytical.fractional.rouse_cv_mid(t, alpha, b, N, D, min_modes=1000)

Velocity autocorrelation of midpoint of a rouse polymer.

Weber Phys Rev E 2010, Eq. 33.

wlcsim.analytical.fractional.rouse_cvv_ep(t, delta, p, alpha, b, N, D)

Term in parenthesis in Lampo, BPJ, 2016 Eq. 10.

wlcsim.analytical.fractional.rouse_nondim_(t, delta, n1, n2, alpha, b, N, D)

uses parameters defined by Lampo et al, BPJ, 2016, eq 12

wlcsim.analytical.fractional.tDeltaN(n1, n2, alpha, b, D)

Lampo et al, BPJ, 2016, eq 11

wlcsim.analytical.fractional.tR(alpha, b, N, D)

Lampo et al, BPJ, 2016, eq 8

wlcsim.analytical.fractional.un_rouse_nondim(tOverDelta, deltaOverTDeltaN, alpha, delN=0.001)

Takes a requested choice of parameters to compare to Tom’s calculated values and generates a full parameter list that satisfies those requirements.

uses approximation defined by Lampo et al, BPJ, 2016, eq 12

In Tom’s paper, delN is non-dimensinoalized away into deltaOverTDeltaN, but that only works in the case where the chain is infinitely long. Otherwise, where exactly we choose n1,n2 will affect the velocity correlation because of the effects of the chain ends. While other choices are truly arbitrary (e.g. D, b, N, etc), delN can be used to specify the size of the segment (relative to the whole chain) used to compute the cross correlation.

wlcsim.analytical.fractional.vc(t, delta, beta)

velocity correlation of locus on rouse polymer. beta = alpha/2.

wlcsim.analytical.fractional.vvc_rescaled_theory(t, delta, beta, A, tDeltaN)

velocity cross correlation of two points on rouse polymer.

wlcsim.analytical.fractional.vvc_unscaled_theory(t, delta, beta, A, tDeltaN)

velocity cross correlation of two points on rouse polymer.

wlcsim.data

wlcsim.input

This module “understands” the input format of wlcsim.exe

class wlcsim.input.InputFormat

An enumeration.

class wlcsim.input.ParsedInput(input_file=None, params=None)

Knows how to handle various input file types used by wlcsim simulator over the years, and transparently converts into new parameter naming conventions.

input = ParsedInput(file_name) print(input.ordered_param_names) # see params in order defined print(input.ordered_param_values) # to see values input.write(outfile_name) # write clone of input file

decide_input_format()

Decide between the two input formats we know of. Not too hard, since one uses Fortran-style comments, which we can look out for, and the other uses bash style comments. Further, the former specifies param names and values on separate lines, while the latter specifies them on the same line.

parse_defines_params_file()

Parse file in the format of src/defines.inc. Each line begins with #define WLC_P__[PARAM_NAME] [_a-zA-Z0-9] where WLC_P__[A-Z]* is the parameter name and the piece after the space is the value of the parameter.

TODO:test

parse_lena_params_file()

Lena-style input files have comment lines starting with a “#”. Any other non-blank lines must be of the form “[identifier][whitespace][value]”, where an identifier is of the form “[_a-zA-Z][_a-zA-Z0-9]*”, and a value can be a boolean, float, int or string. They will always be stored as strings in the params dictionary for downstream parsing as needed.

Identifiers, like fortran variables, are interpreted in a case-insensitive manner by the wlcsim program, and so will be store in all-caps within the ParsedInput to signify this.

parse_original_params_file()

Original-style input files have three garbage lines at the top used for commenting then triplets of lines describing one variable each. The first line of the triplet is a counter, for ease of hardcoding input reader, the second line contains the variable name (which is not used by the input reader) in order to make parsing possible outside of the ahrdcoded wlcsim input reader, and the third line contains the value of the parameter itself. The first two lines of each triplet have a fixed form that we use to extract the record numbers and parameter names, but these forms are not used by wlcsim itself, which ignores these lnies completely. Thus, it is possible that the user specified a particular name for a parameter but that name does not match what wlcsim interpreted it as, since wlcsim simply uses the order of the parameters in this file to determine their identities.

parse_params_file()

Parse and populate ParsedInput’s params, ordered_param_names This parser currently understands three file formats: 1) “ORIGINAL” is the input method understood by Andy’s hardcoded input reader in the original code used the Spakowitz lab was founded. 2) “LENA” is a spec using slightly more general input reader written by Elena Koslover while Andy’s student. 3) “DEFINES” is the format of the src/defines.inc file.

write(output_file_name)

writes out a valid input file for wlcsim.exe given parameters in the format returned by read_file

wlcsim.input.correct_param_name(name)

Takes messy param names from different generations of the simulator and converts them to their newest forms so that simulations from across different years can be tabulated together.

wlcsim.input.correct_param_value(name, value)

Some old param names also have new types. This takes messy param names from different generations of the simulator and converts their names and values to their newest forms so that simulations from across different years can be tabulated together.

wlcsim.FrankElastic

The calculation of Frank elastic constants occures in the function get_Frank_values.

wlcsim.FrankElastic.frank.get_Frank_values(gamma, N, gammaNmax=None, laplace_args={})

Returns dictionary of Frank Elastic Values for specified system.

Nondimensionalized Frank elastc constant K’=K*A/(phi_00*L). Maier saupe parameter ether nondimensionalized by L*A*phi_00 or 2*l_p*A*phi_00 for aLAf and a2lpAf respectively. Correlation functions S are also returned.

Parameters:

• gamma (real) – Field strength in kT’s per Kuhn length of chain
• N (real) – How many Kuhn lengths long the polymers are
• gammaNmax (real) – Cutoff to prevent wasting time.
• laplace_args (dict) – Arguements to pass to inverse laplace calculation.

Note: when N=0.0, gamma is assume to be gamma*N, that is kT’s per chain.

wlcsim.FrankElastic.frank.get_a(gamma, N, gammaNmax=None)

Maier-Saupe parameter need to generate field strength gamma.

Parameters:

• gamma (real) – Field strength in kT’s per Kuhn length of chain
• N (real) – How many Kuhn lengths long the polymers are

Note: when N=0.0, gamma is assume to be gamma*N, that is kT’s per chain.

wlcsim.FrankElastic.frank.get_first_pole(gamma, tol=1e-06)

Search the real axes for polls of G[0][0,0]. Returns highest real value.

Parameters:

• gamma (float) – Field strength.
• tol (float) – absolute tolerance of result

wlcsim.FrankElastic.frank.test_Frank_accuracy(N=10, gamma=1000)

Compair result for different laplace_args.

Looks really good at N=10, gamma=25 To within ~0.3% for N=100, gamma=50. Increasing factor is impoartant. At N=100 gamma=100 the factor=400 level gets off by a few percent. At N=10, gamma=100 it is accurate to may descimal points At N=10, gamma=1000 there is a several persent error

Plotting scripts are found in plot_Frank

wlcsim.FrankElastic.plot_frank.Multi_frank_plot(load=None, logy=False)

Plot data from get_Multi_frank_data

Parameters:

• load (string) – File to load
• logy (bool) – Log y axis

wlcsim.FrankElastic.plot_frank.find_gamma(N, aset=None, gammaN=True, gamma_TOL=1e-06, a_TOL=1e-06)

Frank Elastic data for given N and a

Parameters:

• N (float) – Number of Kuhn length
• aset (itrable) – aLAf or a2lpAf values
• gammaN (bool) – If ture use aLAf, else use a2lpAf, amoung other things

wlcsim.FrankElastic.plot_frank.get_Multi_frank_data(log_gamma=False, npts=48, gammaN=False, maxgamma=35, saveas='Frank.p', nthreads=25, Nset=None)

Calculate Frank elastic constants for a range of gamma and N.

Parameters:

• log_gamma (bool) – Space gamma logrithmeicly
• npts (int) – number gamma values
• N (real) – Kuhn lenths in polymer
• Where to save results. (saveas=) –
• gammaN (bool) – Devide gamma by N. I.e. nondimensionalized by polyer length
• nthreads (int) – Number of threads to use
• Nset (iterable) – Where to save results

wlcsim.FrankElastic.plot_frank.get_Nsweep_data(saveas='Nsweep.p', MultiProcess=True, nthreads=25, Ns=None)

Calculate Frank Elastic Constants for sweep over N

Parameters:

• Where to save results. (saveas=) –
• nthreads (int) – Number of threads to use
• MultiProcess (bool) – Use many threads
• Ns (iterable) – N values

wlcsim.FrankElastic.plot_frank.get_data_for_Splot(name, nthreads=25, N=0.05, npts=45)

Get data for plot_S

Parameters:

• name (string) – Where to save result.
• nthreads (int) – Number of computational threads
• N (real) – Kuhn lenths in polymer
• npts (int) – number gamma values

wlcsim.FrankElastic.plot_frank.multi_process_fun(inputs)

pass through for muiltiproccessing

wlcsim.FrankElastic.plot_frank.plotNsweep(load)

Plot results from get_Nsweep_data

Parameters:load (string) – File to load from.

wlcsim.FrankElastic.plot_frank.plot_Frank(log_gamma=True, MultiProcess=True, npts=24, N=1.0, gammaN=False, maxgamma=35, nthreads=25)

Calculate Frank Elastic data for range of gamma values.

Parameters:

• log_gamma (bool) – Space gamma logrithmeicly
• MultiProcess (bool) – Use multiple threads
• npts (int) – number gamma values
• N (real) – Kuhn lenths in polymer
• gammaN (bool) – Devide gamma by N. I.e. nondimensionalized by polyer length
• nthreads (int) – Number of threads to use

wlcsim.FrankElastic.plot_frank.plot_S(name)

Plot S over diffetn gamma values. :param name: File name to load data from. :type name: string

wlcsim.FrankElastic.plot_frank.seperate_bend_twist_splay(load, gammaN=False, logy=False)

Plot data from get_Multi_frank_data seperatly

Parameters:

• load (string) – File to load
• logy (bool) – Log y axis

There are also a number of helper funtions,

wlcsim.FrankElastic.invlaplace.invLaplace(p_values, G_values, L, reduction=0.0)

Numerical Inverse Laplace from \(G(p) \to G(L)\).

Parameters:

• p_values (array_like) – 1D complex array of p values (i.e. the contour)
• G_values (array_like) – 1D complex array G(p)
• L (float) – Polymer length (In Kuhn lengths)
• reduction (float) – multiply answer by exp(-reduction) for better numerics

wlcsim.FrankElastic.invlaplace.invLaplace_path(path, cutts, G_values, L, reduction=0.0)

Inverse laplace transform based on path.

Parameters:

• path (ndarray) – complex path values
• cutts (dict) – specified by path
• G_values –

wlcsim.FrankElastic.invlaplace.make_path(factor=100, scale=1.0, lambda_offset=0.1, width=2.0, depth=2.0, nwing=100, nside=30, nback=20, maxp=6000, pole_offset=0.0)

Path for complex integral that draws a half rectangle about the origin.

Parameters:

• factor (int) – Number of points to caltulate (sort of)
• scale (float) – Scale of path. Recomended: around 10.0/N
• lambda_offset (float) – distance from complex axis
• width (float) – width of rectangle
• depth (float) – depth of rectangle
• nwing (int) – number of points in wing before factor
• nside (int) – number of points to a side before factor
• back (int) – number of points on the back before factor
• maxp (float) – approximation for infinity
• pole_offset (float) – add to depth

wlcsim.FrankElastic.stonefence.Gmll_matrix

“Matrox of propagators between starting and ending l value.

Parameters:

• Wplus (numpy array) – Result of Wplus_vec for same m, p
• Wminus (numpy array) – Reult of Wminus_vec for same m, p
• Am (numpy array) – Result of Am_vec for same m
• PgammaB (numpy array) – Result of PgammaB_vec for same m, p, gamma
• gamma (float) – alignment strength, in kT’s per Kuhn length
• m (int) – z component of agular momentum quantum number

Returns:

An ORDER_L x ORDER_L numpy matrix with propagators that use Maier-Saupe steps to get from l0 to lf.

wlcsim.FrankElastic.stonefence.PgammaB_vec

P - gamma beta

Returns:vector with index ell

wlcsim.FrankElastic.stonefence.precalculate_data(p, gamma, m_values=[0])

Precalculate W_plus, W_minus, W_pm, and G_m_ll

Parameters:

• p (complex) – laplace conjugate of path length
• gamma (real) – aligning l=2 (Maier-Saupe) field strength
• m_values (list) – list of integer m values to precalculate for

wlcsim.FrankElastic.yyyreal.YR(m, ell, phi, theta)

Real spherical harmonic function.

wlcsim.FrankElastic.yyyreal.YYY(m1, L, M, m2, ORDER_L=50)

Same as YYYreal

wlcsim.FrankElastic.yyyreal.YYY_0mm(m, ORDER_L)

Y_{ell1,0} Y_{2,m} Y_{ell2,m}

wlcsim.FrankElastic.yyyreal.YYY_l1(m1, m, m2, ORDER_L)

Y_{l1,m1} Y_{1,m} Y{l2,m2}

wlcsim.FrankElastic.yyyreal.YYYreal(m1, L, M, m2, ORDER_L)

YYYreal_{l_1,l_2} = int dvec{u} Y_{l_1, m1} Y_{L,M} Y_{l2,m2}

Not able to handle all combiniations of inputs. Can only handle Y_{l_1,0} Y_{2,M} Y_{l_2,M} or Y_{l_1,M} Y_{2,M} Y_{l_2,0} or Y_{l_1,m_1} Y_{1,M} Y{l_2,m_2} or Y_{l_1,m_1} Y_{0,0} Y{l_2,m_2}

wlcsim.FrankElastic.yyyreal.ind_sphere(fun)

Integrate function on sphere.

Parameters:fun (callable) – fun(theta, phi)

wlcsim.FrankElastic.yyyreal.intYYYnum(l1, m1, L, M, l2, m2)

Numerical integral over 3 spherical harmics