-"""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.
-"""
-import numpy as np
-from numba import jit
-# import mpmath
-
-from functools import lru_cache
-from pathlib import Path
-import os
-
-
-[docs]@jit
-def rouse_mode(p, n, N=1):
- """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)"""
- p = np.atleast_1d(p)
- phi = np.sqrt(2)*np.cos(p*np.pi*n/N)
- phi[p == 0] = 1
- return phi
-
-
-[docs]@jit(nopython=True)
-def rouse_mode_coef(p, b, N, kbT=1):
- """k_p: Weber Phys Rev E 2010, after Eq. 18."""
- # alternate: k*pi**2/N * p**2, i.e. k = 3kbT/b**2
- return 3*np.pi**2*kbT/(N*b**2)*p**2
-
-
-[docs]@jit(nopython=True)
-def kp_over_kbt(p : float, b : float, N : float):
- """k_p/(k_B T) : "non-dimensionalized" k_p is all that's needed for most
- formulas, e.g. MSD."""
- return (3*np.pi*np.pi)/(N*b*b) * p*p
-
-
-[docs]@jit(nopython=True)
-def linear_mid_msd(t, b, N, D, num_modes=20000):
- """
- modified from Weber Phys Rev E 2010, Eq. 24.
- """
- msd = np.zeros_like(t)
- for p in range(1, num_modes+1):
- k2p_norm = kp_over_kbt(2*p, b, N)
- msd += (1 / k2p_norm) * (1 - np.exp(-k2p_norm * (D / N) * t))
- return 12 * msd + 6 * D * t / N
-
-
-[docs]@jit(nopython=True)
-def gaussian_G(r, N, b):
- """Green's function of a Gaussian chain at N Kuhn lengths of separation,
- given a Kuhn length of b"""
- r2 = np.power(r, 2)
- return np.power(3/(2*np.pi*b*b*N), 3/2)*np.exp(-(3/2)*r2/(N*b*b))
-
-
-[docs]@jit(nopython=True)
-def gaussian_Ploop(a, N, b):
- """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"""
- Nb2 = N*b*b
- return spycial.erf(a*np.sqrt(3/2/Nb2)) - a*np.sqrt(6/np.pi/Nb2)/np.exp(3*a*a/2/Nb2)
-
-
-@jit(nopython=True)
-def _cart_to_sph(x, y, z):
- r = np.sqrt(x*x + y*y + z*z)
- if r == 0.0:
- return 0.0, 0.0, 0.0
- phi = np.arctan2(y, x)
- theta = np.arccos(z/r)
- return r, phi, theta
-
-
-def confined_G(r, rp, N, b, a, n_max=100, l_max=50):
- # first precompute the zeros of the spherical bessel functions, since our
- # routine to do so is pretty slow
- if l_max >= 86: # for l >= 86, |m| >= 85 returns NaN from sph_harm
- raise ValueError("l_max > 85 causes NaN's from scipy.special.sph_harm")
- if confined_G.zl_n is None or n_max > confined_G.zl_n.shape[1] \
- or l_max > confined_G.zl_n.shape[0]:
- confined_G.zl_n = spherical_jn_zeros(l_max, n_max)
- # some aliases
- spherical_jn = scipy.special.spherical_jn
- Ylm = scipy.special.sph_harm
- zl_n = confined_G.zl_n[:l_max+1,:n_max]
- # convert to spherical coordinates
- r = np.array(r)
- if r.ndim == 1:
- x, y, z = r
- xp, yp, zp = rp
- # elif r.ndim == 2:
- # x = r[:,0]
- # y = r[:,1]
- # z = r[:,2]
- # xp = rp[:,0]
- # yp = rp[:,1]
- # zp = rp[:,2]
- else:
- raise ValueError("Don't understand your R vectors")
- r, phi, theta = _cart_to_sph(x, y, z)
- rp, phip, thetap = _cart_to_sph(xp, yp, zp)
- # compute terms based on indexing (ij in meshgrid to match zl_n)
- l, n = np.meshgrid(np.arange(l_max+1), np.arange(n_max), indexing='ij')
- ln_term = 2*np.exp(-b**2/6 * (zl_n/a)**2 * N)
- ln_term = ln_term/(a**3 * spherical_jn(l+1, zl_n)**2)
- ln_term = ln_term*spherical_jn(l, zl_n/a*r)*spherical_jn(l, zl_n/a*rp)
- # l,m terms
- m = np.arange(-l_max, l_max+1)
- l, m = np.meshgrid(np.arange(l_max+1), np.arange(-l_max, l_max+1))
- lm_term = Ylm(m, l, phi, theta)*Ylm(m, l, phip, thetap)
- lm_mask = np.abs(m) <= l
- lm_term[~lm_mask] = 0
- # now broadcast and sum
- G = np.sum(ln_term[None,:,:] * lm_term[:,:,None])
- return G
-confined_G.zl_n = None
-
-
-[docs]def linear_mscd(t, D, Ndel, N, b=1, num_modes=20000):
- r"""
- 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
- :math:`k_BT/\xi` for :math:`\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 : (M,) np.array<float>
- result
- """
- mscd = np.zeros_like(t)
-
- k1 = 3 * np.pi ** 2 / (N * (b ** 2))
- sum_coeff = 48 / k1
- exp_coeff = k1 * D / N
- sin_coeff = np.pi * Ndel / N
-
- for p in range(1, num_modes+1, 2):
- mscd += (1/p**2) * (1 - np.exp(-exp_coeff * (p ** 2) * t)) \
- * np.sin(sin_coeff*p)**2
-
- return sum_coeff * mscd
-
-
-[docs]def ring_mscd(t, D, Ndel, N, b=1, num_modes=20000):
- r"""
- 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
- :math:`k_BT/\xi` for :math:`\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 : (M,) np.array<float>
- result
- """
- mscd = np.zeros_like(t)
-
- k1 = 12 * np.pi ** 2 / (N * (b ** 2))
- sum_coeff = 48 / k1
- exp_coeff = k1 * D / N
- sin_coeff = 2 * np.pi * Ndel / N
-
- for p in range(1, num_modes+1):
- mscd += (1 / p ** 2) * (1 - np.exp(-exp_coeff * p ** 2 * t)) \
- * np.sin(sin_coeff * p) ** 2
- return sum_coeff * mscd
-
-
-[docs]def end_to_end_corr(t, D, N, num_modes=10000):
- """Doi and Edwards, Eq. 4.35"""
- mscd = np.zeros_like(t)
- tau1 = N**2/(3*np.pi*np.pi*D)
- for p in range(1, num_modes+1, 2):
- mscd += 8/p/p/np.pi/np.pi * np.exp(-t*p**2 / tau1)
- return N*mscd
-
-import numpy as np
-
-
-[docs]def r2_ave(length_kuhn, dimensions=3):
- r"""
- r2_ave - Calculate the average end-to-end distance squared :math:`\langle R^{2} \rangle / (2 l_{p})^{2}`
- for the wormlike chain model
-
- Parameters
- ----------
- length_kuhn : float (array)
- The length of the chain in Kuhn lengths
- dimensions : int
- The number of dimensions (default to 3 dimensions)
-
- Returns
- -------
- r2 : float (array)
- The mean-square end-to-end distance for the wormlike chain model (non-dimensionalized by :math:`2 l_{p})`
-
- Notes
- -----
- See Mehraeen, et al, Phys. Rev. E, 77, 061803 (2008). (Ref [Mehraeen2008]_)
- """
- r2 = 2 * (length_kuhn / (dimensions - 1)
- - (1 - np.exp(-(dimensions - 1) * length_kuhn)) / (dimensions - 1) ** 2)
-
- return r2
-
-
-[docs]def rg2_ave(length_kuhn, dimensions=3):
- r"""
- rg2_ave - Calculate the radius of gyration
- :math:`\langle \vec{R}_{G}^{2} \rangle / (2 l_{p})^{2}` for the wormlike chain model
-
- Parameters
- ----------
- length_kuhn : float (array)
- The length of the chain in Kuhn lengths
- dimensions : int
- The number of dimensions (default to 3 dimensions)
-
- Returns
- -------
- rg2 : float (array)
- The mean-square radius of gyration for the wormlike chain model (non-dimensionalized by :math:`2 l_{p})`
-
- Notes
- -----
- See Mehraeen, et al, Phys. Rev. E, 77, 061803 (2008). (Ref [Mehraeen2008]_)
- """
- rg2 = 2 * (length_kuhn / (6 * (dimensions - 1)) - 1 / (2 * (dimensions - 1) ** 2)
- + length_kuhn ** -1 / (dimensions - 1) ** 3
- - length_kuhn ** -2 * (1 - np.exp(-(dimensions - 1) * length_kuhn)) / (dimensions - 1) ** 4)
-
- return rg2
-
-
-[docs]def rz4_ave(length_kuhn, dimensions=3):
- r"""
- rz4_ave - Calculate the 4th moment of the end-to-end distribution
- :math:`\langle R_{z}^{4} \rangle / (2 l_{p})^{4}` for the wormlike chain model
-
- Parameters
- ----------
- length_kuhn : float (array)
- The length of the chain in Kuhn lengths
- dimensions : int
- The number of dimensions (default to 3 dimensions)
-
- Returns
- -------
- rz4 : float (array)
- The mean-square end-to-end distance for the wormlike chain model (non-dimensionalized by :math:`2 l_{p})`
-
- Notes
- -----
- See Mehraeen, et al, Phys. Rev. E, 77, 061803 (2008). (Ref [Mehraeen2008]_)
- """
-
- # Calculate the coefficients from the ladder operations on the hyperspherical harmonics (dimensions)
- a1 = np.sqrt(1 / dimensions)
- a2 = np.sqrt(2 * (dimensions - 1) / ((dimensions + 2) * dimensions))
-
- # Calculate the contributions from the 2 contributing stone-fence diagrams
- diagram1 = (length_kuhn ** 2 / (2 * (dimensions - 1) ** 2) - 2 * length_kuhn / (dimensions - 1) ** 3
- + 3 / (dimensions - 1) ** 4
- - length_kuhn * np.exp(-(dimensions - 1) * length_kuhn) / (dimensions - 1) ** 3
- - 3 * np.exp(-(dimensions - 1) * length_kuhn) / (dimensions - 1) ** 4)
- diagram1 *= a1 ** 4
-
- diagram2 = (length_kuhn / (2 * dimensions * (dimensions - 1) ** 2)
- - (5 * dimensions - 1) / (4 * dimensions ** 2 * (dimensions - 1) ** 3)
- + np.exp(-(dimensions - 1) * length_kuhn) * (
- length_kuhn / ((dimensions + 1) * (dimensions - 1) ** 2)
- + (dimensions + 3) / ((dimensions - 1) ** 3 * (dimensions + 1) ** 2))
- + np.exp(-2 * dimensions * length_kuhn) / (4 * dimensions ** 2 * (dimensions + 1) ** 2))
- diagram2 *= a1 ** 2 * a2 ** 2
-
- rz4 = 24 * (diagram1 + diagram2)
-
- return rz4
-
-
-from wlcstat.util.wlc_poles_residues import *
-import numpy as np
-import scipy.special as sp
-
-
-[docs]def eval_poles_and_residues(k_val, mu, mu_zero_only=True, lam_zero_only=True, dimensions=3):
- r"""
- eval_poles_and_residues - Evaluate the poles and the residues for a given value of the
- Fourier vector magnitude :math:`K`
-
- Parameters
- ----------
- k_val : float
- The value of the Fourier vector magnitude :math:`K`
- mu : int
- Value of the mu parameter (:math:`z`-component of the angular momentum)
- mu_zero_only : boolean
- Determines whether poles and residues are determined for non-zero :math:`\mu` (default True)
- lam_zero_only : boolean
- Determines whether residues are determined for non_zero :math:`\lambda` (default True)
- dimensions : int
- The number of dimensions (default to 3 dimensions)
-
- Returns
- -------
- poles : complex float
- Evaluated poles for the given :math:`K` and :math:`\mu`
- residues : complex float
- Evaluated residues for the given :math:`K` and :math:`\mu`
-
- Notes
- -----
- See Mehraeen, et al, Phys. Rev. E, 77, 061803 (2008). (Ref [Mehraeen2008]_)
- """
-
- poles = eval_poles(k_val, mu, dimensions)
- residues = eval_residues(k_val, mu, poles, lam_zero_only, dimensions)
-
- return poles, residues
-
-
-[docs]def gwlc_r(r_val, length_kuhn, dimensions=3, alpha_max=25, k_val_max=1e5, delta_k_val_max=0.1):
- r"""
- gwlc_r - Evaluate the orientation-independent Green's function for the wormlike chain model
-
- Parameters
- ----------
- r_val : float (array)
- The values of the end-to-end distance :math:`r = R/L` to be evaluated
- length_kuhn : float (array)
- The length of the chain in Kuhn lengths
- dimensions : int
- The number of dimensions (default to 3 dimensions)
- alpha_max : int
- Maximum number of poles evaluated (default 25)
- k_val_max : float
- Cutoff value of :math:`K` for numerical integration
- delta_k_val_max : float
- Maximum value of the integration step size
-
- Returns
- -------
- gwlc : float
- The orientation-independent Green's function [size len(r_val) x len(length_kuhn)]
-
- Notes
- -----
- See Mehraeen, et al, Phys. Rev. E, 77, 061803 (2008). (Ref [Mehraeen2008]_)
- """
-
- delta_k_val = min(delta_k_val_max, 2 * np.pi / np.max(length_kuhn) / 10)
-
- # Eliminate 0 and 1 from the r_val array
- r_val[r_val == 0] = 1e-10
- r_val[r_val == 1] = 1-1e-10
-
- # Initialize the Green's function
- if type(length_kuhn) == float or type(length_kuhn) == int:
- gwlc = np.zeros((len(r_val)), dtype=type(1+1j))
- else:
- gwlc = np.zeros((len(r_val), len(length_kuhn)), dtype=type(1+1j))
-
- tolerance = 1e-15
-
- k_val_output = k_val_max / 100
- k_val = delta_k_val
- int_count = 0
- contains_nan = False
- while k_val <= k_val_max and not contains_nan:
- int_count += 1
-
- poles = eval_poles(k_val, 0, dimensions, alpha_max)
- residues = eval_residues(k_val, 0, poles, True, dimensions, alpha_max, alpha_max)
-
- if int_count == 1:
- int_coef = 55 / 24
- elif int_count == 2:
- int_coef = -1 / 6
- elif int_count == 3:
- int_coef = 11 / 8
- else:
- int_coef = 1
-
- for alpha in range(0, alpha_max + 1):
- gkwlc_kval = residues[alpha] * np.exp(poles[alpha] * length_kuhn)
- if type(length_kuhn) == float or type(length_kuhn) == int:
- integrand = (int_coef * k_val ** (dimensions / 2)
- * sp.jv(dimensions / 2 - 1, k_val * r_val * length_kuhn) * gkwlc_kval)
- else:
- integrand = (int_coef * k_val ** (dimensions / 2)
- * sp.jv(dimensions / 2 - 1, k_val * np.outer(r_val, length_kuhn))
- * np.outer(np.ones((len(r_val)), dtype=type(1+1j)), gkwlc_kval))
- if not contains_nan:
- contains_nan = np.isnan(integrand).any()
- if not contains_nan:
- gwlc += integrand
- else:
- print("Encountered NaN at k_val = " + str(k_val))
-
- k_val += delta_k_val
- if k_val >= k_val_output:
- k_val_output += k_val_max / 100
- print("Current k_val = " + str(k_val) + " with delta_k_val = " + str(delta_k_val)
- + " and k_val_max = " + str(k_val_max))
- print(np.max(abs(np.exp(poles * np.min(length_kuhn)))))
-
- if np.min(abs(np.exp(poles * np.min(length_kuhn)))) < tolerance and alpha_max > 10 and k_val > 2840:
- alpha_max -= 1
- print("Reducing alpha_max to " + str(alpha_max) + " at k_val = "
- + str(k_val) + " with delta_k_val = " + str(delta_k_val))
- print(np.max(abs(np.exp(poles * np.min(length_kuhn)))))
-
- if np.max(abs(np.exp(poles * np.min(length_kuhn)))) < tolerance:
- print("Achieved accuracy of " + str(tolerance) + " at k_val = "
- + str(k_val) + " with delta_k_val = " + str(delta_k_val))
- print(np.max(abs(np.exp(poles * np.min(length_kuhn)))))
- k_val = k_val_max + delta_k_val
-
- gwlc *= delta_k_val / (2 * np.pi) ** (dimensions / 2) * np.outer(r_val ** (-dimensions / 2 + 1),
- length_kuhn ** (dimensions / 2 + 1))
-
- return gwlc
-
-
-from wlcstat.util.wlc_poles_residues import *
-from wlcstat.util.wlc_vertex import *
-import numpy as np
-
-
-[docs]def eval_structure_factor(k_val_vector, length_kuhn, dimensions=3, alpha_max=25):
- r"""
- eval_structure factor - Evaluate the structure factor for the wormlike chain model
-
- Parameters
- ----------
- k_val_vector : float (array)
- The value of the Fourier vector magnitude :math:`K`
- length_kuhn : float (array)
- The length of the chain in Kuhn lengths
- dimensions : int
- The number of dimensions (default to 3 dimensions)
- alpha_max : int
- Maximum number of poles evaluated (default 25)
-
- Returns
- -------
- structure_factor : float (vector)
- Structure factor for the wormlike chain model for every k_val in k_val_vector
-
- Notes
- -----
- See Mehraeen, et al, Phys. Rev. E, 77, 061803 (2008). (Ref [Mehraeen2008]_)
- """
- if type(length_kuhn) == float or type(length_kuhn) == int:
- structure_factor = np.zeros((len(k_val_vector)), dtype=type(1+1j))
- else:
- structure_factor = np.zeros((len(k_val_vector), len(length_kuhn)), dtype=type(1+1j))
-
- for ind_k_val in range(0, len(k_val_vector)):
- k_val = k_val_vector[ind_k_val]
-
- poles = eval_poles(k_val, 0, dimensions, alpha_max)
- residues = eval_residues(k_val, 0, poles, True, dimensions, alpha_max)
- residue_zero, ddp_residue_zero = eval_residue_zero(k_val, dimensions)
-
- if type(length_kuhn) == float or type(length_kuhn) == int:
- structure_factor[ind_k_val] = (length_kuhn * residue_zero + ddp_residue_zero)
- for alpha in range(0, alpha_max):
- structure_factor[ind_k_val] += (np.exp(poles[alpha] * length_kuhn) *
- residues[alpha] / poles[alpha] ** 2)
- structure_factor[ind_k_val] *= 2 / length_kuhn ** 2
- else:
- for ind_length in range(0, len(length_kuhn)):
- structure_factor[ind_k_val, ind_length] = (length_kuhn[ind_length] * residue_zero + ddp_residue_zero)
-
- for alpha in range(0, alpha_max):
- structure_factor[ind_k_val, ind_length] += (np.exp(poles[alpha] * length_kuhn[ind_length]) *
- residues[alpha] / poles[alpha] ** 2)
-
- structure_factor[ind_k_val, ind_length] *= 2 / length_kuhn[ind_length] ** 2
-
- return structure_factor
-
-
-