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/(N1) 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
[source]¶ 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
[source]¶ Unfortunately, by “cleaning up” this function, we also make it 2x slower.
#worthit

wlcsim.bd.rouse.
jit_srk1
[source]¶ Faster version of wlcsim.bd.rouse.rouse using jit.
N=101,L=100,b=1,D=1
takes about 3.5min to run whent=np.linspace(0, 1e5, 1e7+1)
addingt_save
does not slow function downOur 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 fastk to rouselike 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')[source]¶ Get the fullpolymer 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 crossover 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 threeregime MSD behavior can be very easily visualized in loglog 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 loglog space), have slopes 1, 1/2, and 1, respectively, and their “offsets” (yintercept 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 shorttime behavior to get \(D\) either, even supposing you manage to measure short enough times to see the crossover 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 3d 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 yintercept of the MSD in loglog 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 3d 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 simulationcentric 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 loglog space intersects \(6 \frac{\hat{D}}{\tilde{N}} t\) at \(t = t_R\), you get the above form for the subRouse 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)[source]¶ 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>)[source]¶ 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 singlelocus 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 freedraining 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(i1, t)) + R(t)\]where \(\xi = k_BT/D\), \(k = 3k_BT/b^2\), \(b\) is the Kuhn length of the polymer, \(D\) is the selfdiffusion coefficient of a bead, and \(R(t)/\xi\) is a deltacorrelated stationary Gaussian process with mean zero and \(\langle R(t) R(t')/\xi^2 \rangle = 2DI\delta(tt')\).
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 homologouslylinked 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
[source]¶ “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)[source]¶ 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 interbead segments), and one tenth of the beads bound to each other, we might get the following Poissonspace 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 101105 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 106117 are linked sequentially.
Similarly, beads 18 and 118 (and 27 and 125) are linked, with beads 118125 being linked sequentially.

wlcsim.bd.homolog.
rouse
(N, FP, L, b, D, Aex, rx, ry, rz, t, t_save=None, tether_list=None)[source]¶ BD simulation of two homologous yeast chromosomes in meiosis.
An arbitrary ellipticalshaped 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()
(andf_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 shorttime regular diffusion or Rouselike 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)[source]¶ 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.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 justintime compileable when possible.

wlcsim.bd.forces.
f_elas_homolog_rouse
[source]¶ 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
[source]¶ Compute spring forces on single, linear rouse polymer.
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
[source]¶ center bead by default
for 1e4long time arrays, this takes ~1030s on my laptop

wlcsim.bd.runge_kutta.
ou
(x0, t, k_over_xi, D, method=<function rk4_thermal_lena>)[source]¶ 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)[source]¶ 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)[source]¶ 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)[source]¶ 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, stretchableshearable 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)[source]¶ 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)[source]¶ TODO: test

wlcsim.bd.tsswlc.
sswlc
(N, L, lp, t, t_save)[source]¶ Simulate a stretchableshearable 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)[source]¶ 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 endtoend 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 endtoend distance, the function will ValueError.

wlcsim.analytical.rouse.
end2end_distance_gauss
(r, b, N, L)[source]¶ 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)[source]¶ Doi and Edwards, Eq. 4.35

wlcsim.analytical.rouse.
gaussian_G
[source]¶ Green’s function of a Gaussian chain at N Kuhn lengths of separation, given a Kuhn length of b

wlcsim.analytical.rouse.
gaussian_Ploop
[source]¶ 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
[source]¶ “nondimensionalized” k_p is all that’s needed for most formulas, e.g. MSD.
Type: k_p/(k_B T)

wlcsim.analytical.rouse.
linear_mscd
(t, D, Ndel, N, b=1, num_modes=10000)[source]¶ 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)[source]¶ 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)[source]¶ Cvv^delta(t) for infinite polymer.
Lampo, BPJ, 2016 Eq. 16.

wlcsim.analytical.rouse.
rouse_mode
[source]¶ 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.
test_rouse_msd_line_approx
()[source]¶ 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 thatmsd_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, somsd_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.01321183e07
for N = 1e8+1; L = 174; b = 150; D = 166;, we have msd_approx(t)  msd(t) = 5.2889657(3)e05
for N = 1e8+1; L = 174; b = 15; D = 166;, we have msd_approx(t)  msd(t) = 5.2889657(3)e06
for N = 1e8+1; L = 17.4; b = 15; D = 166;, we have msd_approx(t)  msd(t) = 5.2889657(3)e07
for N = 1e8+1; L = 17.4; b = 15; D = 16.6;, we have msd_approx(t)  msd(t) = 5.2889657(3)e07 (no change)
for N = 1e7+1; L = 17.4; b = 15; D = 16.6;, we have msd_approx(t)  msd(t) = 5.2889657(3)e06 (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)[source]¶ Velocity autocorrelation of a fractionallydiffusing particle. Weber, Phys Rev E, 2010 (Eq 32)

wlcsim.analytical.fractional.
frac_discrete_cv
(t, delta, alpha, D)[source]¶ Discrete velocity autocorrelation of a fractionallydiffusing particle. Weber, Phys Rev E, 2010 (Eq 33)

wlcsim.analytical.fractional.
frac_discrete_cv_normalized
(t, delta, alpha)[source]¶ Normalized discrete velocity autocorrelation of a fractionallydiffusing 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)[source]¶ 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=1e05, atol=1e08, force_convergence=True)[source]¶ Velocity crosscorrelation 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)[source]¶ Weber Phys Rev E 2010, Eq. 24.

wlcsim.analytical.fractional.
frac_rouse_mode_corr
(p, t, alpha, b, N, D)[source]¶ Weber Phys Rev E 2010, Eq. 21.

wlcsim.analytical.fractional.
rouse_cv_mid
(t, alpha, b, N, D, min_modes=1000)[source]¶ 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)[source]¶ Term in parenthesis in Lampo, BPJ, 2016 Eq. 10.

wlcsim.analytical.fractional.
rouse_nondim_
(t, delta, n1, n2, alpha, b, N, D)[source]¶ uses parameters defined by Lampo et al, BPJ, 2016, eq 12

wlcsim.analytical.fractional.
un_rouse_nondim
(tOverDelta, deltaOverTDeltaN, alpha, delN=0.001)[source]¶ 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 nondimensinoalized 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)[source]¶ velocity correlation of locus on rouse polymer. beta = alpha/2.
wlcsim.data¶
wlcsim.input¶
This module “understands” the input format of wlcsim.exe

class
wlcsim.input.
ParsedInput
(input_file=None, params=None)[source]¶ 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
()[source]¶ Decide between the two input formats we know of. Not too hard, since one uses Fortranstyle 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
()[source]¶ Parse file in the format of src/defines.inc. Each line begins with #define WLC_P__[PARAM_NAME] [_azAZ09] where WLC_P__[AZ]* is the parameter name and the piece after the space is the value of the parameter.
TODO:test

parse_lena_params_file
()[source]¶ Lenastyle input files have comment lines starting with a “#”. Any other nonblank lines must be of the form “[identifier][whitespace][value]”, where an identifier is of the form “[_azAZ][_azAZ09]*”, 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 caseinsensitive manner by the wlcsim program, and so will be store in allcaps within the ParsedInput to signify this.

parse_original_params_file
()[source]¶ Originalstyle 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
()[source]¶ 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.


wlcsim.input.
correct_param_name
(name)[source]¶ 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)[source]¶ 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={})[source]¶ 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)[source]¶ MaierSaupe 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=1e06)[source]¶ 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)[source]¶ 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)[source]¶ 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=1e06, a_TOL=1e06)[source]¶ 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)[source]¶ 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)[source]¶ 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)[source]¶ 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)[source]¶ pass through for muiltiproccessing

wlcsim.FrankElastic.plot_frank.
plotNsweep
(load)[source]¶ 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)[source]¶ 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)[source]¶ 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)[source]¶ 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)[source]¶ 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)[source]¶ 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)[source]¶ 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
[source]¶ “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 MaierSaupe steps to get from l0 to lf.

wlcsim.FrankElastic.stonefence.
precalculate_data
(p, gamma, m_values=[0])[source]¶ Precalculate W_plus, W_minus, W_pm, and G_m_ll
Parameters:  p (complex) – laplace conjugate of path length
 gamma (real) – aligning l=2 (MaierSaupe) field strength
 m_values (list) – list of integer m values to precalculate for

wlcsim.FrankElastic.yyyreal.
YYYreal
(m1, L, M, m2, ORDER_L)[source]¶ 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}