# -*- coding: utf-8 -*-
# written by Ralf Biehl at the Forschungszentrum Jülich ,
# Jülich Center for Neutron Science (JCNS-1)
# Jscatter is a program to read, analyse and plot data
# Copyright (C) 2015-2026 Ralf Biehl
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import inspect
import os
import numbers
import numpy as np
from scipy import linalg as la
import scipy.constants
from .. import dataArray as dA
from .. import dataList as dL
from .. import formel
try:
from ..libs import fscatter
useFortran = True
except ImportError:
fscatter = None
useFortran = False
__all__ = ['solveOptimizedRouseZimm', 'linearChainORZ','multiArmStarORZ', 'ringChainORZ']
pi = np.pi
_path_ = os.path.realpath(os.path.dirname(__file__))
#: Planck constant in µeV*ns
h = scipy.constants.Planck / scipy.constants.e * 1E15 # µeV*ns
#: h/2π reduced Planck constant in µeV*ns
hbar = h/2/pi # µeV*ns
kb = 1.3806503e-23 # J/K in SI units
[docs]
def solveOptimizedRouseZimm(A, reducedfriction=0.25):
r"""
Solve optimized Rouse-Zimm (ORZ) approximation to the diffusion equation of a polymer in solution.
From [1]_ :
A generalization of some work by Bixon on the theoretical foundations of the Rouse-Zimm model in
polymer solution dynamics. In particular, a procedure is described for constructing the "best possible"
Rouse-Zimm model for an arbitrary polymer, starting from the equilibrium distribution of polymer conformations
and using either Kirkwood's generalized diffusion dynamics or stochastic dynamics.
The method is based on an application of linear response theory to the calculation
of certain time correlation functions for polymer dynamics.
Parameters
----------
A : [N,N] array
Structural matrix.
reducedfriction : float
Reduced friction :math:`\zeta_r= \zeta/6\pi\eta_s l` in matrix of the preaveraged hydrodynamic interactions
:math:`H_{ij}=\delta_{ij} + \zeta_r \langle l/R_{ij}\rangle (1-\delta_{ij})`.
The default 0.25 is in accordance with experimental data for theta solvents [2]_.
- =0 : free draining limit, no HI, Rouse dynamics. H is identity matrix that :math:`HA=A`.
- >0 : with HI typically <0.5. For h>0.5 we may get negative eigenvalues.
[3]_ : The value of :math:`\zeta_r = 0.25` ensures that the matrix [H] is
positive definite and does not show any instabilities related to
the appearance of negative unphysical eigenvalues. Such
eigenvalues occur when the parameter :math:`\zeta_r` exceeds a critical
value :math:`\zeta_r^*` which corresponds to the non-overlapping
condition :math:`r/b=\frac{bead radius}{bond length} ≤ 0.5` for the monomers.
[6]_ : The abrupt change in :math:`\phi_{max}` at r/b~O.43 indicates that a large perturbation takes
place in the modes. Therefore, this value has to be considered an upper limit to the strength of
the hydrodynamic interaction. More restrictive conditions, physically reasonable though arbitrary,
fall in the range of r/b lower than 0.43 but do not apparently give rise to anomalies in the modes.
In the application to intrinsic viscosity in ideal solvents, presented in Sec. III,
we shall therefore take r/b in the range 0-0.40.
Returns
-------
evalHA, evecHA, A, loverR, Rg2_red, Rij2_red : array's
evalHA : 1D array
Eigenvalues :math:`\lambda_a` of :math:`HA`.
evecHA : 2D array
Eigenvectors :math:`Q_{ia}` of :math:`HA`.
mu : 1D array
Diagonal elements :math:`\mu_a = (Q^TAQ)_{aa}`
loverR : 2D array
Adimensional mean inverse distance matrix equ. 25 in [4]_ :
:math:`\langle l/R_{ij}\rangle =l(6/\pi)^{1/2}(\langle R^2_{ij}\rangle )^{-1/2}` .
Rg2_red : float
Reduced radius of gyration (equ 3.2 in [2]_)
:math:`R_g^2/l^2 = \sum_{ij} \langle R^2_{ij}\rangle / (2N^2)`
Rij2_red : 2D array
Reduced second moments for the distance between any two chain atoms equ. 3.1 in [2]_ :
.. math:: \langle R^2_{ij}\rangle /l^2= \sum_{k=1}^{N} (Q_{ik} - Q_{jk})^2/\lambda_k
A : 2D array
Dimesionless structural matrix :math:`A = M[1:,:].T * U * M[1:,:]`
(or force constant matrix as :math:`3kT/l^2 A` or connectivity matrix)
Notes
-----
In the **ORZ approximation** to the Kirkwood configurational diffusion equation,
the bead coordinates of a chain of N beads (or monomers) of friction coefficient :math:`\zeta` jointed by equal
links of mean square length :math:`l^2` and forceconstant :math:`\kappa` obey a Langevin equation of the form
.. math:: \zeta \frac{\partial}{\partial t} R_i(t) + \kappa \sum_{j} (HA)_{ij} R_j(t) = v_i(t) \\
\frac{\partial}{\partial t} R_i(t) + \sigma \sum_{j} (HA)_{ij} R_j(t) = v_i(t)
were :math:`v_i` is the random velocity, :math:`\sigma=\kappa/\zeta=3kT/l^2\zeta` as bond rate constant [4]_
and bead coordinates :math:`R_i`.
:math:`\zeta=6\pi\eta_0l` is the bead friction with the surrounding/solvent.
The structural matrix :math:`A` depends on the actual shape of the polymer (e.g. linear, star or ring).
The matrix H is the preaveraged hydrodynamic matrix with elements
.. math:: H_{ij} = \delta_{ij} + \zeta_r\langle l/R_{ij}\rangle (1-\delta_{ij})
with :math:`R_{ij}` as distance between beads i and j, reduced friction per chain atom
:math:`\zeta_r=\zeta/6\pi\eta_0l` and :math:`\eta_0` as solvent viscosity.
The structural matrix A of a **linear** polymer is
.. math:: A = M^T \left( \begin{array}{cc} 0 & 0 \\ 0 & U \end{array} \right) M = M[1:,:]^T U M[1:,:]
with the transfer matrix M in dimsionless form (the first row is not needed)
.. math:: M = \left ( \begin{array}{cccc} 1/N & 1/N & ... & 1/N \\
-1 & 1 & 0 & 0..0 \\
0 & -1 & 1 & 0..0 \\
... \end{array} \right )
In the first dimension is the center of mass that is basically not needed.
Elements of :math:`M_{ij}` are −1 if the bond vector :math:`l_i` starts at monomer `i` and
+1 if bond vector :math:`l_i` points to bead `i`, else 0.
U can be retrieved from (depends on the model: freely jointed chain (FJC), free rotating chain (FRC), RIS,...)
.. math:: U^{-1}_{ij} = \langle l_i\cdot l_j\rangle /l^2
- Freely jointed linear chain and bead spring model:
.. math:: \langle l_i\cdot l_j\rangle /l^2 = \delta_{ij} \Rightarrow \langle R^2_{ij}\rangle = l^2 |i-j|
- Free rotating linear chain with bond angles :math:`g=-cos(\theta)` (:math:`\theta=\pi` is rigid rod)
.. math:: \langle l_i\cdot l_j\rangle/l^2 = g^{|i-j|} \Rightarrow
\langle R^2_{ij}\rangle = l^2 |i-j| [\frac{1+g}{1-g} -\frac{2g}{|i-j|}\frac{1-g^{|i-j|}}{(1-g)^2}]
Persistence length is :math:`nl=l/(1-g)` and :math:`\langle R^2_{ij}\rangle` from [5]_.
For :math:`g=0` we yield the FJC.
The **solution of a ORZ model** can be described by normal modes and corresponding eigenvalues
:math:`Q_a, \lambda_a` of the matrix :math:`HA` yielding normal coordinates :math:`\Xi_k`,
mode relaxation times :math:`\tau_a`
and mean square displacements :math:`\langle \xi^2_a \rangle` and
.. math:: \xi_k(t) = Q_{ki} R_i(t)
.. math:: \tau_a=(\sigma \lambda_a)^{-1}
.. math:: \langle \xi^2_a\rangle = l^2 \mu_a^{-1}
where :math:`\mu_a= (Q^TAQ)_{aa}`. For free draining (:math:`\zeta_r=0`) :math:`\mu_a=\lambda_a`
According to the usual Zimm notation, the mode :math:`a=0` describes the translational mode of
the center of resistance, always characterized by a constant eigenvector :math:`Q_{k,0} = N^{-1/2}` and
a zero relaxation rate :math:`\lambda_0 = 0`
The dynamic correlation between beads i and j is [1]_ :
.. math:: \langle |R_i(t) -R_j(0)|^2\rangle &= l^2 \sum_{a=1}^{N-1} \mu_a^{-1} [|Q_{ia}|^2+|Q_{ja}|^2 -
(Q_{ia}Q_{ka}^* + Q_{ia}Q_{ka}^*) exp(-\sigma\lambda_at)] \\
&= l^2 \sum_{a=1}^N \mu_a^{-1} [|Q_{ia}|^2+|Q_{ja}|^2 - 2Q_{ia}Q_{ka} exp(-\sigma\lambda_at)] \\
&= l^2 d_{ij}^2
where the first refers to complex eigenvectors [2]_ and the later to real eigenvectors [4]_.
The dynamic structure factor (measured by DLS or NSE) is [4]_
(here extended by bead scattering amplitudes :math:`f_{i}`):
.. math:: S(q,t)/S(q,0) &= \frac{1}{F(q)} \sum_{ij}^{N} f_{i}f_{j} exp(-q^2/6\langle |R_i(t) -R_j(0)|^2\rangle ) \\
&= \frac{1}{F(q)} exp(-q^2Dt) \sum_{ij}^{N} f_{i}f_{j }exp(-q^2l^2/6 \cdot d_{ij}^2)
with the form factor :math:`F(q)` and (equ. 8 in [4]_)
.. math:: d_{ij}^2 = \sum_{a=1}^N \mu_a^{-1} [|Q_{ia}|^2+|Q_{ja}|^2 - 2Q_{ia}Q_{ka} exp(-\sigma\lambda_at))]
.. note::
One should recognize the similarity in the description to normal modes of biomacromolecules reprectively proteins
as described in the Ornstein-Uhlenbeck process :py:func:`jscatter.bio.intScatFuncOU`.
The difference is the explicit usage of coordinates in :py:func:`jscatter.bio.intScatFuncOU`
while in the ORZ model statistical averages of :math:`d_{ij}` without explicit coordinates are used.
We yield the static form factor as for :math:`t \rightarrow 0`
.. math:: F(Q) &= S(q,0) =& \sum_{ij}^N f_{i}f_{j}exp(-\frac{q^2}{6}\langle R_{ij}^2\rangle ) \\
\hat{F}(Q) &= & \frac{1}{N^2}\sum_{ij}^N exp(-\frac{q^2}{6}\langle R_{ij}^2\rangle )
where the second line is the normalized form factor equ. 2 in [4] with :math:`f_i=1`.
The translational diffusion coefficient is equ. 24 in [4]:
.. math:: D = kT/N\zeta [1 + \frac{\zeta_r}{N}\sum_{ij} (1-\delta_{ij}) \langle l/R_{ij}\rangle ]
The first cumulant observed in an experiment at short times
(initial slope :math:`\Omega(q)=-(d/dt)ln(S(q,t))|_{t=0}=\Omega(q)`) is equ. 13 in [4]_ for :math:`f_i=1`:
.. math:: \Omega(q) = q^2D_{cum} = q^2 \frac{\sigma l^2}{3N} \frac{1}{\hat{F}(q)} [1 + \frac{\zeta_r}{N}
\sum_{ij} (1-\delta_{ij}) \langle l/R_{ij}\rangle exp(-q^2l^2/6 \cdot d_{ij}^2(0) )]
The prefactor is :math:`\sigma l^2/(3N)=kT/N\zeta` as in the equation above.
Taking into account the scattering amplitudes :math:`f_i` we yield with the NOT normalized form factor :math:`F(q)`
.. math:: \Omega(q) = q^2 D_{cum} = q^2 \frac{\sigma l^2}{3} \frac{1}{F(q)F(0)}
[\sum_{i=j}f_i^2 + \sum_{i\neq j}
f_i f_j \zeta_r \langle l/R_{ij}\rangle exp(-q^2l^2/6 \cdot d_{ij}^2(0) )]
Here we see that specific bead/arm contributions can be suppressed if the respective beads
are matched to the surrounding solvent.
Nevertheless, the overall tranlational diffusion will not change by matching.
The dynamic intrinsic viscosity can be calculated from the relaxation times equ. 21 in [5]_ :
.. math:: [\eta(\omega)] = \frac{N_akT}{M\eta_0} \sum_{a} \frac{\tau_a}{1+i\omega\tau_a}
with M as molecular weight and :math:`N_a` as Avogadro constant.
.. note::
It should be pointed out that the matrix :math:`U^{-1}` geht highly singular if the model
is applied to more rigid chains or rodlike parts resulting in negative eigenvalues.
**Internal friction**
Internal friction :math:`\zeta_{int}` between neigboring beads can be included in analogy to Soranno [7]_
.. math:: \zeta \frac{\partial}{\partial t} R_i(t) + \kappa \sum_{j} (HA)_{ij} R_j(t) +
\zeta_{int} \frac{\partial}{\partial t} (HA)_{ij} R_j(t) = v_i(t)
with the additional friction term. Using the above eigenvectors of :math:`HA` this leads to
.. math:: \zeta \frac{\partial}{\partial t} \xi_k(t) + \kappa \lambda_k \xi_k(t) +
\zeta_{int} \frac{\partial}{\partial t} \lambda_k \xi_k(t) &= w_k(t) \\
[(\zeta +\zeta_{int} \lambda_k)\frac{\partial}{\partial t}\xi_k(t) = - \kappa \lambda_k \xi_k(t) + w_k(t)
and result in same eigenvectors but changed relaxation times
:math:`\tau_{a,int}=(\sigma \lambda_a)^{-1} + \tau_{int} = \tau_a + \tau_{int}`
with :math:`\tau_{int} = \zeta_{int}/\kappa` as presented by Sorrano for standard Rouse/Zimm model [7]_.
The additional internal friction is here independent of the modes and has a stronger effect on higher modes.
In above corrrelation functions we need to change
.. math:: &exp(-\sigma\lambda_a t)=exp(- \frac{t}{\tau_a}) \rightarrow \\
&exp(-\frac{\sigma \lambda_a}{1+\sigma \lambda_a \tau_{int}}t)=exp(-\frac{t}{\tau_a+\tau_{int}})
Here :math:`\tau_a` contains implicitly the mode dependence :math:`\tau_{zp}=\tau_zp^{-3\nu}` for Zimm or
:math:`\tau_{rp}=\tau_zp^{-2}` for Rouse like systems (see :py:func:`finiteRouse`, :py:func:`finiteZimm`)
but the additional mode independent internal friction is the same.
Correspondingly in the cumulant a correction is needed
:math:`\langle l/R_{ij}\rangle \rightarrow \langle l/R_{ij}\rangle \sum_a 1/(1+\sigma\lambda_a \tau_{int})`
Examples
--------
Example to reproduce Fig 6 of [2]_ (only left subplot) but in simpler model (not RIS)
just using bead-spring (FJC) and free rotating chain models for linear polymers:
Non-vanishing bond correlations increase size.
Draining increases relaxation times but not the msd of the movements.
::
import jscatter as js
import numpy as np
p = js.grace(1.9,1.)
p.multi(1,3,hgap=0.25)
p[0].yaxis(label=r'\xl\f{}\sa')
p[1].yaxis(label=r'\xt\f{}\sa\N\xs',scale='log',ticklabel=['power',0])
p[2].yaxis(label=r'\x<z\f{}\S2\N\sa\N>/l\S2',scale='log',ticklabel=['power',0])
p[0].xaxis(label=r'a')
p[1].xaxis(label=r'a',scale='log')
p[2].xaxis(label=r'a',scale='log')
# matrices for linear chain
N=100
M = np.diag([1.]*(N+1)) + np.diag([-1.]*N,-1)
U = np.diag([1.]*N)
A = M[1:,:].T @ U @ M[1:,:]
ev, evec, mu = js.dynamic.solveOptimizedRouseZimm(A, reducedfriction=0.25)[:3]
a = np.r_[1:ev.shape[0]]
p[0].plot(a, ev[1:],le='FJC partial draining')
p[1].plot(a, 1/ev[1:],le='FJC partial draining')
p[2].plot(a, 1/mu[1:],le='FJC partial draining')
ev, evec, mu = js.dynamic.solveOptimizedRouseZimm(A, reducedfriction=0)[:3]
a = np.r_[1:ev.shape[0]]
p[0].plot(a, ev[1:],le='FJC free draining')
p[1].plot(a, 1/ev[1:],le='FJC free draining')
p[2].plot(a, 1/mu[1:],le='FJC free draining')
# matrices for linear chain but non-vanishing bond correlation
costheta = 0.65
i,j = np.indices([N,N])
Uinv = costheta**np.abs(i-j)
U = np.linalg.inv(Uinv)
A = M[1:,:].T @ U @ M[1:,:]
ev, evec, mu = js.dynamic.solveOptimizedRouseZimm(A, reducedfriction=0.25)[:3]
a = np.r_[1:ev.shape[0]]
p[0].plot(a, ev[1:],le='FRC partial draining')
p[1].plot(a, 1/ev[1:],le='FRC partial draining')
p[2].plot(a, 1/mu[1:],le='FRC partial draining')
ev, evec, mu = js.dynamic.solveOptimizedRouseZimm(A, reducedfriction=0)[:3]
a = np.r_[1:ev.shape[0]]
p[0].plot(a, ev[1:],le='FRC free draining')
p[1].plot(a, 1/ev[1:],le='FRC free draining')
p[2].plot(a, 1/mu[1:],le='FRC free draining')
p[0].subtitle('eigenvalue spectra')
p[1].subtitle('relaxation times')
p[2].subtitle('mean square displacements')
p[1].title('model FJC + FRC free draining and partial draining')
p[1].legend(x=1.4,y=1,charsize=0.8)
# p.save(js.examples.imagepath+'/ORZeigenvalue.jpg',size=(1.9,1.),dpi=200)
.. image:: ../../examples/images/ORZeigenvalue.jpg
:align: center
:width: 80 %
:alt: ORZ eigenvalue and more
Correctness of the solution can be verified by comparing to Table I of [6]_.
The parameter `h` in [6]_ corresponds to reducedfriction :math:`=r/b= (\pi/(3N))^{1/2}h \approx 0.1023 h`
Today the eigenvalue problem can be solved directly and is more accurate.
::
import jscatter as js
import numpy as np
# matrices for linear chain
N=100
M = np.diag([1.]*(N+1)) + np.diag([-1.]*N,-1)
U = np.diag([1.]*N)
A = M[1:,:].T @ U @ M[1:,:]
# for h=1 or 2 compare to :math:`\lambda_k` (second and 4th column) and :math:`\mu_k` for h=0 to exact solution.
ev, evec, mu = js.dynamic.solveOptimizedRouseZimm(A, reducedfriction=0.2)[:3]
ev[-1] # ~ 2.66
mu[-1] # ~ 3.99
ev, evec, mu = js.dynamic.solveOptimizedRouseZimm(A, reducedfriction=0.1)[:3]
ev[-1] # ~ 3.33
mu[-1] # ~ 3.99
References
----------
.. [1] Theoretical basis for the Rouse‐Zimm model in polymer solution dynamics.
Zwanzig, R. The Journal of Chemical Physics 60, 2717–2720 (1974) https://doi.org/10.1063/1.1681433
.. [2] A hierarchy of models for the dynamics of polymer chains in dilute solution.
Perico, A., Ganazzoli, F. & Allegra, G.
The Journal of Chemical Physics 87, 3677–3686 (1987). https://doi.org/10.1063/1.452966
.. [3] Intramolecular relaxation dynamics in semiflexible dendrimers.
Kumar, A. & Biswas, P. Journal of Chemical Physics 134, (2011). https://doi.org/10.1063/1.3598336
.. [4] Static and Dynamic Structure Factors for Star Polymers in θ Conditions.
Guenza, M. & Perico, A.
Macromolecules 26, 4196–4202 (1993). https://doi.org/10.1021/ma00068a020
.. [5] Optimized Rouse–Zimm theory for stiff polymers
M. Bixon; R. Zwanzig
J. Chem. Phys. 68, 1896–1902 (1978) https://doi.org/10.1063/1.435916
.. [6] Dynamics of chain molecules. I. Solutions to the hydrodynamic equation and intrinsic viscosity.
Perico, A., Piaggio, P. & Cuniberti, C.
The Journal of Chemical Physics 62, 4911–4918 (1975). https://doi.org/10.1063/1.430404
.. [7] Internal friction in an intrinsically disordered protein - Comparing Rouse-like models with experiments
A. Soranno, F. Zosel, H. Hofmann
J. Chem. Phys. 148, 123326 (2018) http://aip.scitation.org/doi/10.1063/1.5009286
"""
N = A.shape[0]
# eigenvalues of static problem (H=Identity) give Rij2
# [2] equ 3.1 and below
# also Kumar et al J. Chem. Phys. 138, 104902 (2013) equ 3 But not mentioning using static problem
# numpy : The column evec[:, i] is the normalized eigenvector to the eigenvalue eval[i]. evals are sorted
evalA, evecA= la.eigh(A)
Rij2_red = fscatter.dynamic.eigvector2rij2(evalA, evecA) # not including l**2
# [2] equ 3.2 + 3.3 without lo**2
# Rg2_red = Rg2/ l**2 = np.sum(Rij2_red)/N**2/2 or
Rg2_red = np.sum(1/evalA[1:]) / N
# get loverR, catch diagonal zeros
np.fill_diagonal(Rij2_red, 1)
loverR = (6 / np.pi)**0.5 * Rij2_red ** -0.5
np.fill_diagonal(Rij2_red, 0)
np.fill_diagonal(loverR, 0)
if reducedfriction > 0:
# H = hydrodynamic matrix, [1] equ 9
H = np.diag([1]*N) + reducedfriction * loverR
# print(f'H has {np.sum(la.eigvals(H) <0)} neg eigenvalues. NOT positive definite ')
# if A is symmetric and H symmetric and positive definite then H@A can be diagonalized
evalHA, evecHA = la.eigh(H @ A)
# calc mode amplitude factor (evec.T @ A @ evec) is diagonal up to numerical precision
muHA = np.diag(evecHA.T @ A @ evecHA)
return evalHA, evecHA, muHA, loverR, Rg2_red, Rij2_red, A
else:
# Rouse = no HI interaction => H ==identity
# calc mode amplitude factor (evec.T @ A @ evec) is diagonal up to numerical precision
muA = np.diag(evecA.T @ A @ evecA)
return evalA, evecA, muA, loverR, Rg2_red, Rij2_red, A
@formel.memoize(maxsize=4)
def _linearStructuralMatrix(N, costheta):
"""
Create for a star the bond correlation matrix U^-1 and bead-to-bond vector transformation matrix M
This part is in reduced units that bond length is not used.
"""
# see references in main function
# bond correlation matrix U_ij^-1 = <l_i*l_j>/l**2 ; with N-1 bonds
if isinstance(costheta, numbers.Number):
if costheta==0:
# identity for uncorrelated bonds , freely jointed chain no correlation except self correlation
Uinv = np.diag([1]*(N-1))
U = Uinv # np.linalg.inv(Uinv)
else:
# free rotating chain FRC with <li*lj> /l**2 = cos(theta)^(|j-i|) or Prod[cos(theta_i)]
# eg. https://dasher.wustl.edu/bio5357/readings/rubinstein-chapter2.pdf 2.20
# [1] equ. 30-33 ff
i,j = np.indices([N-1,N-1])
# same arm => difference in indices gives distance
Uinv = costheta**np.abs(i-j) # |i-j| < n_arm
Uinv[Uinv < 0.001] = 0
U = np.linalg.inv(Uinv)
elif isinstance(costheta, (list, np.ndarray)):
assert len(costheta) == N-2, 'costheta should be list of len N-2.'
assert np.all(np.array(costheta) < 1) & np.all(
0 <= np.array(costheta)), 'costheta should be 0 <= cos(theta) < 1 !'
i,j = np.indices([N-1,N-1])
Uinv = fscatter.dynamic.cumcos(costheta, i, j)
Uinv[Uinv < 0.001] = 0
U = np.linalg.inv(Uinv)
else:
# for later cos theta depends on bead
raise TypeError('costheta should be float or list of len N-1.')
# bead-to-bond vector transformation matrix [1] equ 25 particular for star
# indices in paper run from 1 to n , here 0 to N-1
M = np.diag([0.]+[1.]*(N-1)) + np.diag([-1.]*(N-1),-1)
# [2] equ 2.7 A.T@U @ A with a = M[1:,:]
A = M[1:,:].T @ U @ M[1:,:]
return A
[docs]
def linearChainORZ(t, q, N=100, l=1, fa=None, Dcm=None, Dcmfkt=None,
viscosity=1, costheta=0., T=293, reducedfriction=0.25, tintern=0):
r"""
Dynamics of a linear polymer chain using optimized Rouse-Zimm approximation (ORZ).
The linear dynamic structure factor is calculated in analogy to the star described by Guenza [1]_.
The linear chain dynamics is decribed within the more general optimized Rouse-Zimm (ORZ) approximation
introduced by Bixon and Zwanzig [3]_ [4]_ . See :py:func:`solveOptimizedRouseZimm`.
We extend this here using bead scattering lengths for partial matching and allow individual bond correlations.
To speedup fitting use the memoize function as described in :py:func:`~.formel.memoize`.
Parameters
----------
t : array
timepoints in units ns.
q : array
Scattering vectors in units 1/nm.
N : int
Number of beads .
l : float, default = 0.38 nm (amino acid)
Bond length or Kuhn length in units nm.
fa : None, list of float
Scattering length of bead/monomer :math:`fa_i`. Can be used to match parts of the star to the solvent.
- None : Equal 1 for all beads.
- list: length N as scattering length of beads in sequence
Dcm : float
Center of mass diffusion in units nm**2/ns.
If `None` the calculated value D_ORZ is used.
Dcmfkt : array 2xN, function
Function f(q) or array with [qi, f(qi) ] as correction for Dcm like Diff = Dcm*f(q).
e.g. for inclusion of structure factor or hydrodynamic function with f(q)=H(Q)/S(q).
Missing values are interpolated.
costheta : float, list of float, 0 <= costheta < 1
Cos of bond correlation angle :math:`\langle \vec{l_i} \cdot \vec{l_j} \rangle /l^2 = cos(theta)`
between normalized bonds :math:`l_i` with :math:`0 \le cos(\theta) \le 1`.
- costheta = 0 : FJC (freely jointed chain) model, no bond correlation, Rouse dynamics, No HI .
- float :math:`0 < cos(\theta) \le 1` : FRC (free rotating chain). With
:math:`R_{ee} = Nl^2C_{\infty}=Nl^2 \frac{1+cos(\theta)}{1-cos(\theta)}` .
- list of float length N-2 : FRC with individual :math:`cos(\theta_i)` for each pair of N-2 bonds.
- ``costheta=([0.1]*24+[0.7]*24)`` for stretched beginning and flexible end of 50 beads.
- ``costheta=0.8 * np.cos(np.pi*np.r_[:1:(N-2)*1j])**2 + 0.1`` for flexible center and stretched ends.
For fitting encapsulate the model in a function where you parametrize your model for costheta .
- During fits use ``limits(costheta=[None,None,0.001,0.999])`` to avoid singular matrices.
e.g. (from [5]_ p. 53) :math:`C_{\infty} | Kuhn length [A] | cos(theta)`
- polyisoprene 4.6 | 8.2 | 0.783
- polyethylene oxide 6.7 | 11 | 0.85
- polyethylene 7.4 | 14 | 0.865
- atactic polystyrene 9.5 | 18 | 0.895
viscosity : float, default=1 (H2O@20C )
:math:`\eta` in units cPoise=mPa*s e.g. water :math:`visc(T=293 K) =1 mPas`
T : float, default 273+20
Temperature in Kelvin.
reducedfriction : float, default 0.25
Reduced friction :math:`\zeta_r = \zeta/6\pi\eta_sl` in hydrodynamic tensor H.
- =0 : free draining limit, no HI, Rouse dynamics.
- >0 : with HI typically <0.5. For h>0.5 we get negative eigenvalues.
See :py:func:`solveOptimizedRouseZimm`.
During fits use ``limits(reducedfriction=[0,0.43,0.,0.5])`` to avoid singular matrices.
tintern : float>0, default 0
Relaxation time due to internal friction between neighboring beads in ns.
Returns
-------
sqt : dataList
Intermediate scattering function of a star for given q.
- [times; Sqt; Sqt only diffusion; Sqt cumulant]
- columnname = 't;Sqt;Sqt_inf;Sqt_cum'
- .q : scattering vector
- .D_ORZcum : diffusion coefficinet in initial slope (cumulant)
- .D_ORZ : translational diffusion ORZ model
- .Dcm : used center of friction/mass diffusion
- .Fq : form factor
- .Fq_inf : form factor t=inf
- .beadfriction : bead friction
- .bondrateconstant : bondrateconstant in 1/ns
- .moderelaxationtimes : used moderelaxationtimes in ns [4] equ 2.17
- .mode_rmsd : used mode rmsd as math:`l/\mu^{0.5} # [4] equ 2.19
- .reducedfriction : reducedfriction
- .costheta : costheta
- .eigenvalues : all evals
- .mu : all mu
- .l : bondlength l
- .N : Number of beads
- .Rg : radius of gyration in nm
- .Rg_red : reduced radius of gyration
Notes
-----
See :py:func:`solveOptimizedRouseZimm` for a description of the ORZ model with respective parameters and
the dynamic structure factor :math:`S(q,t)/S(q,0)` .
Here we use a linear chain with N beads and set elements :math:`U^{-1}_{ij} <0.001 \rightarrow 0`.
The inverse of the static bond correlation matrix :math:`U_{ij}^{-1} = \langle l_i\cdot l_j\rangle /l^2`
in dimesionless form is
.. math:: U_{ij}^{-1} &= \delta_{i,j} &\text{ for uncorrelated bonds } \\
&= \prod_{n=i}^{j} g_n \; &\text{ for individual } g_i \text{ including constant g}
The transfer matrix M is
.. math:: M = \delta_{i,j} - \delta_{i+1,j}
Examples
--------
Here we examine how changig stiffness influences dynamics.
::
import jscatter as js
import numpy as np
q= np.r_[0.01,0.1:2:0.2]
t = np.r_[0:1:0.02,1:20:1,20:100:5]
def stiffendschain(t, q, N, l=0.5, cosmin =0.05,cosmax=0.8,rf=0.25):
costheta = (cosmax-cosmin) * np.cos(np.pi*np.r_[:1:(N-2)*1j])**2 + cosmin
sqt = js.dynamic.linearChainORZ(t, q, N, l=0.5, costheta = costheta , reducedfriction=rf)
return sqt
p = js.grace()
p.xaxis(label='t / ns')
sqt = stiffendschain(t, q, 100, l=1, cosmin =0.1,cosmax=0.5)
for c,sq in enumerate(sqt,1):
p.plot(sq.X,sq._Sqt,li=[1,1,c],sy=0,le=f'{sq.q:.2f} nm\\S-1')
sqt = stiffendschain(t, q, 100, l=1, cosmin =0.1,cosmax=0.1)
for c,sq in enumerate(sqt,1):
p.plot(sq.X,sq._Sqt,li=[3,1,c],sy=0,le='')
p.yaxis(label='S(Q,t),S(Q,0)',scale='log',min=0.01,max=1)
p.legend(charsize=0.7)
p.title('chain with stiff ends ')
p.subtitle('solid: stiff ends; broken: flexible ends')
# p.save(js.examples.imagepath+'/ORZ_linearstiffends.png',size=(1.5,1.5),dpi=300)
.. image:: ../../examples/images/ORZ_linearstiffends.png
:align: center
:width: 60 %
:alt: ORZ eigenvalue and more
References
----------
.. [1] Static and Dynamic Structure Factors for Star Polymers in θ Conditions.
Guenza, M. & Perico, A.
Macromolecules 26, 4196–4202 (1993). https://doi.org/10.1021/ma00068a020
.. [2] A Local Approach to the Dynamics of Star Polymers.
Guenza, M., Mormino, M. & Perico, A.
Macromolecules 24, 6168–6174 (1991). https://doi.org/10.1021/ma00023a018
.. [3] Theoretical basis for the Rouse‐Zimm model in polymer solution dynamics.
Zwanzig, R. The Journal of Chemical Physics 60, 2717–2720 (1974) https://doi.org/10.1063/1.1681433
.. [4] A hierarchy of models for the dynamics of polymer chains in dilute solution.
Perico, A., Ganazzoli, F. & Allegra, G.
The Journal of Chemical Physics 87, 3677–3686 (1987). https://doi.org/10.1063/1.452966
.. [5] Rubinstein, M. & Colby, R. H. Polymer Physics. (OUP Oxford, 2003).
.. [6] Intramolecular relaxation dynamics in semiflexible dendrimers.
Kumar, A. & Biswas, P. Journal of Chemical Physics 134, (2011). https://doi.org/10.1063/1.3598336
"""
if fa is None:
fa = np.ones(N)
assert len(fa) == N, 'fa should be of length N.'
fa = np.array(fa)
q =np.atleast_1d(q)
t0 = np.r_[0,t]
# bead friction coefficient
friction = 6 * np.pi * (viscosity*1e-3) * (l * 1e-9) # l in nm; viscosity in Pa*s
# bond rate constant brc, sometimes named W
brc = (3 * kb * T / (l*1e-9)**2 / friction) * 1e-9 # l in nm # brc is 1/s => *1e-9 in 1/ns
reducedfriction = max(min(reducedfriction, 1), 0)
# create correction for diffusion
if Dcmfkt is not None:
if formel._getFuncCode(Dcmfkt):
# is already an interpolation function
Dcmfunktion = Dcmfkt
elif np.shape(Dcmfkt)[0] == 2:
Dcmfunktion = lambda qq: dA(Dcmfkt).interp(qq)
else:
raise TypeError('Shape of Dcmfkt is not 2xN!')
else:
# by default no correction
Dcmfunktion = lambda qq: 1.
# create bond matrix and transfer matrix for a star
A = _linearStructuralMatrix(N, costheta)
# compute Eigenvectors and Eigenvalues for ORZ in reduced units
evals, evec, mu, loverR, Rg2_red, Rij2_red2, _ = solveOptimizedRouseZimm(A, reducedfriction)
if np.any(evals[1:] < 0):
raise UserWarning(f'There are {np.sum(evals[1:]<0)} negative eigenvalues in ORZ solution. '
f'reducedfriction should be smaller.')
# >95% of computing time in this call for N=100 that already uses omp
# array dim [Q values x time values+2]
# first is sqt(0), last element is sum in equation 13 for initial cumulant, second last t=inf
# first eval is COM diffusion and large eval and large mu dont contribute much
lowev = np.r_[[False], (brc * evals[1:] < 15) | (1/mu[1:] > max(0.001/mu[1:]))]
sqt = fscatter.dynamic.sqtnonlinearpolymer(evec[:,1:], evals[1:], mu[1:], fa, q, t0, l, brc, loverR, tintern)
# [1] equ 24 for center of mass diffusion
Dcm_ORZ = kb * T / N / friction * 1e9 # in nm**2/ns
if reducedfriction>0:
Dcm_ORZ *= (1 + reducedfriction/N * np.sum(loverR))
results = dL()
for i, sqti in enumerate(sqt):
# calc initial slope diffusion as first cumulant/q**2 [1] equ 13
D_cum_ORZ = brc * l**2/3/sqti[0] * (np.sum(fa**2) + reducedfriction * sqti[-1])
if Dcm is None:
Dcm = Dcm_ORZ
results.append(np.c_[t0[1:],
np.exp(-q[i]**2 * Dcm * Dcmfunktion(q[i]) * t0[1:]) * sqti[1:-2]/sqti[0],
np.exp(-q[i]**2 * Dcm * Dcmfunktion(q[i]) * t0[1:]) * sqti[-2]/sqti[0],
np.exp(-q[i]**2 * (D_cum_ORZ - Dcm + Dcm * Dcmfunktion(q[i])) * t0[1:]) ].T)
results[-1].setColumnIndex(iey=None)
results[-1].columnname = 't;Sqt;Sqt_inf;Sqt_cum'
results[-1].viscosity = viscosity
results[-1].q = q[i]
results[-1].D_ORZcum = D_cum_ORZ
results[-1].D_ORZ = Dcm_ORZ
results[-1].Dcm = Dcm
results[-1].Fq = sqti[0] / np.sum(fa)**2 # form factor t=0
results[-1].Fq_inf = sqti[-2] / np.sum(fa)**2 # form factor t=inf
results[-1].beadfriction = friction
results[-1].bondrateconstant = brc # in 1/ns
results[-1].moderelaxationtimes = 1 / (brc * evals[lowev]) + tintern # [4] equ 2.17
results[-1].mode_rmsd = l / mu[lowev]**0.5 # [4] equ 2.19
results[-1].reducedfriction = reducedfriction
results[-1].costheta = costheta
results[-1].eigenvalues = evals
results[-1].mu = mu
results[-1].l = l
results[-1].tintern = tintern
results[-1].N = evec.shape[0]
results[-1].Rg = l * Rg2_red**0.5
results[-1].Rg_red = Rg2_red**0.5
results[-1].modelname = inspect.currentframe().f_code.co_name
if len(results) == 1:
return results[0]
return results
@formel.memoize(maxsize=4)
def _starStructuralMatrix(f_arm, n_arm, costheta):
"""
Create for a star the bond correlation matrix U^-1 and bead-to-bond vector transformation matrix M
This part is in reduced units that bond length is not used.
"""
# see references in main function
# f_arms number of arms
# n_arms number of beads in arm
N = f_arm * n_arm + 1 # total number of beads, first index is center of the star
# bond correlation matrix U_ij^-1 = <l_i*l_j>/l**2 ; with N-1 bonds
if isinstance(costheta, numbers.Number):
if costheta==0:
# identity for uncorrelated bonds , freely jointed chain no correlation except self correlation
Uinv = np.diag([1]*(N-1))
U = Uinv # np.linalg.inv(Uinv)
else:
# free rotating chain FRC with <li*lj> /l**2 = cos(theta)^(|j-i|) or Prod[cos(theta_i)]
# eg. https://dasher.wustl.edu/bio5357/readings/rubinstein-chapter2.pdf 2.20
# [1] equ. 30-33 ff
i,j = np.indices([N-1,N-1])
Uinv = np.zeros([N-1,N-1])
# same arm => difference in indices gives distance
same = (i//n_arm == j // n_arm)
Uinv[same] = costheta**np.abs(i-j)[same] # |i-j| < n_arm
# different arms => each distance to center + center arm correlation
# similar to equ [1] 32,33 but a=1/(1-f) and indices i,j start at 0 not 1 like paper [1]
different = (i//n_arm != j // n_arm)
if f_arm > 2:
# center is symmetric
Uinv[different] = (1/(f_arm -1) * costheta**(i%n_arm + j%n_arm ))[different]
else:
# for linear case there is no special about the center, 1/(f-1) makes a singular matrix
Uinv[different] = (costheta * costheta**(i % n_arm + j % n_arm))[different]
Uinv[Uinv < 0.001] = 0
U = np.linalg.inv(Uinv)
elif isinstance(costheta, (list, np.ndarray)):
#raise UserWarning('not yet implemented')
assert len(costheta) == n_arm-1, 'costheta should be list of len n_arm-1.'
assert np.all(np.array(costheta) < 1) & np.all(
0 <= np.array(costheta)), 'costheta should be 0 <= cos(theta) < 1 !'
cumcos = fscatter.dynamic.cumcos
i,j = np.indices([N-1,N-1])
Uinv = np.zeros([N-1,N-1])
same = (i//n_arm == j // n_arm)
Uinv[same] = cumcos(costheta, i%n_arm, j%n_arm)[same] # calcs for to much, cut these by [same]
# different arms => each distance to center + center arm correlation
# similar to equ [1] 32,33 but a=1/(1-f) and indices i,j start at 0 not 1 like paper [1]
different = (i//n_arm != j // n_arm)
z0 = np.zeros_like(i)
if f_arm > 2:
# center is symmetric, we use 1/(f-1) (see paper [2])
Uinv[different] = (1/(f_arm -1) * cumcos(costheta,z0, j%n_arm) *
cumcos(costheta,z0, i%n_arm))[different]
else:
# for linear case there is no special about the center, 1/(f-1) would make a singular matrix
Uinv[different] = (costheta[0] * cumcos(costheta,z0, j%n_arm) *
cumcos(costheta,z0, i%n_arm))[different]
Uinv[Uinv < 0.001] = 0
U = np.linalg.inv(Uinv)
else:
# for later cos theta depends on bead
raise TypeError('costheta should be float or list of len n_arm.')
# bead-to-bond vector transformation matrix [1] equ 25 particular for star
# indices in paper run from 1 to n , here 0 to N-1
M = np.diag([0.]+[1.]*(N-1))
# M[0,:] = 1/N # this is never accessed or used
for j in np.r_[:f_arm]:
for i in np.r_[1:n_arm]:
ij = n_arm * j + i
M[ij+1,ij] = -1.
for j in np.r_[:f_arm]:
i = 1 + j * n_arm
M[j*n_arm+1,0] = -1.
# [2] equ 2.7 A.T@U @ A with a = M[1:,:]
A = M[1:,:].T @ U @ M[1:,:]
return A
[docs]
def multiArmStarORZ(t, q, f_arm=4, n_arm=10, l=1, fa=None, Dcm=None, Dcmfkt=None,
viscosity=1, costheta=0., T=293, reducedfriction=0.25, tintern=0):
r"""
Dynamics of a symmetric multi arm star of polymer chains using optimized Rouse-Zimm approximation (ORZ).
The star dynamic structure factor is explicitly described by Guenza [1]_ and extended to partially stretched stars
in [2]_.
The star polymer dynamics is decribed within the more general optimized Rouse-Zimm (ORZ) approximation
introduced by Bixon and Zwanzig [3]_ [4]_ . See :py:func:`solveOptimizedRouseZimm`.
We extend this here using bead scattering lengths for partial matching and allow individual bond correlations.
To speedup fitting use the memoize function as described in :py:func:`~.formel.memoize`.
Parameters
----------
t : array
timepoints in units ns.
q : array
Scattering vectors in units 1/nm.
f_arm : int
Number of arms :math:`f_{arm}`.
n_arm : int
Number of beads per arm :math:`n_{arm}` excluding the center common for all arms.
l : float, default = 0.38 nm (amino acid)
Bond length or Kuhn length in units nm.
fa : None, list of float
Scattering length of bead/monomer :math:`fa_i`. Can be used to match parts of the star to the solvent.
- None : Equal 1 for all beads.
- list: length (1 + f_arm * n_arm) as scattering length of beads in sequence
[center, 1..n_arm first arm,1..n_arm second arm,....]
e.g. ``fa=[0] + ([0]*5+[1]*6)*4``
for a 4 arm star of 11 beads per arm with the center and 5 innermost beads matched.
Dcm : float
Center of mass diffusion in units nm**2/ns.
If `None` the calculated value D_ORZ is used.
Dcmfkt : array 2xN, function
Function f(q) or array with [qi, f(qi) ] as correction for Dcm like Diff = Dcm*f(q).
e.g. for inclusion of structure factor or hydrodynamic function with f(q)=H(Q)/S(q).
Missing values are interpolated.
costheta : float, list of float, 0 <= costheta < 0.1
Cos of bond correlation angle :math:`\langle \vec{l_i} \cdot \vec{l_j} \rangle /l^2 = cos(theta)`
between normalized bonds :math:`l_i` with :math:`0 \le cos(\theta) \le 1`.
- costheta = 0 : FJC (freely jointed chain) model, no bond correlation, Rouse dynamics, No HI .
- float :math:`0 < cos(\theta) \le 1` : FRC (free rotating chain). With
:math:`R_{ee} = Nl^2C_{\infty}=Nl^2 \frac{1+cos(\theta)}{1-cos(\theta)}` .
- list of float length N-2 : FRC with individual :math:`cos(\theta_i)` for each pair of N-2 bonds.
- ``costheta=([0.1]*24+[0.7]*24)`` for stretched beginning and flexible end of 50 beads.
- ``costheta=0.8 * np.cos(np.pi*np.r_[:1:(N-2)*1j])**2 + 0.1`` for flexible center and stretched ends.
For fitting encapsulate the model in a function where you parametrize your model for costheta .
- During fits use ``limits(costheta=[None,None,0.001,0.999])`` to avoid singular matrices.
e.g. (from [5]_ p. 53) :math:`C_{\infty} | Kuhn length [A] | cos(theta)`
- polyisoprene 4.6 | 8.2 | 0.783
- polyethylene oxide 6.7 | 11 | 0.85
- polyethylene 7.4 | 14 | 0.865
- atactic polystyrene 9.5 | 18 | 0.895
viscosity : float, default=1 (H2O@20C )
:math:`\eta` in units cPoise=mPa*s e.g. water :math:`visc(T=293 K) =1 mPas`
T : float, default 273+20
Temperature in Kelvin.
reducedfriction : float, default 0.25
Reduced friction :math:`\zeta_r = \zeta/6\pi\eta_sl` in hydrodynamic tensor H.
- =0 : free draining limit, no HI, Rouse dynamics.
- >0 : with HI typically <0.5. For h>0.5 we get negative eigenvalues.
See :py:func:`solveOptimizedRouseZimm`.
During fits use ``limits(reducedfriction=[0,0.43,0.,0.5])`` to avoid singular matrices.
tintern : float>0, default 0
Relaxation time due to internal friction between neighboring beads in ns.
Returns
-------
sqt : dataList
Intermediate scattering function of a star for given q.
- [times; Sqt; Sqt only diffusion; Sqt cumulant]
- columnname = 't;Sqt;Sqt_inf;Sqt_cum'
- .q : scattering vector
- .D_ORZcum : diffusion coefficinet in initial slope (cumulant)
- .D_ORZ : translational diffusion ORZ model
- .Dcm : used center of friction/mass diffusion
- .Fq : normalized form factor
- .Fq_inf : normalized form factor t=inf
- .beadfriction : bead friction
- .bondrateconstant : bondrateconstant in 1/ns
- .moderelaxationtimes : used moderelaxationtimes in ns [4] equ 2.17
- .mode_rmsd : used mode rmsd as math:`l/\mu^{0.5} # [4] equ 2.19
- .reducedfriction : reducedfriction
- .costheta : costheta
- .eigenvalues : all evals
- .mu : all mu
- .l : bondlength l
- .N : Number of beads
- .Rg : radius of gyration in nm
- .Rg_red : reduced radius of gyration
Notes
-----
See :py:func:`solveOptimizedRouseZimm` for a description of the ORZ model with respective parameters and
the dynamic structure factor :math:`S(q,t)/S(q,0)` .
Here we use a symetric star with `f_arm` arms of each `n_arm` beads and a connecting bead
as described by Guenza [1]_ and set elements :math:`U^{-1}_{ij} <0.001 \rightarrow 0`.
The inverse of the static bond correlation matrix :math:`U_{ij}^{-1} = \langle l_i\cdot l_j\rangle /l^2`
in dimesionless form is
.. math:: U_{ij}^{-1} &= \delta_{i,j} &\text{ for uncorrelated bonds } \\
For individual :math:`g_i` including that g all are the same and core :math:`g_0` (indexing start at 1)
.. math::
U_{ij}^{-1} &= \prod_{n=i}^{j} g_n &\text{ for i,j on the same arm} \\
&= g_0 \prod_{n=1}^{j} g_n \prod_{m=1}^{i} g_m &\text{ for i,j on different arms}
The transfer matrix M is (ignoring the not used :math:`M_{1i}`) accordig to [2]_ :
.. math::
\begin{align}
&M_{1,i} &= 1/N &\text{ for }i=1..n_{arm} \\
&M_{i,i} &= 1 &\text{ for }i=2..n_{arm} \\
&M_{i+1,i} &= -1 &\text{ for }i=2..n_{arm},n_{arm}+2..2n_{arm},...,(f_{arm}-1)N+2...f_{arm}n_{arm} \\
&M_{i,1} &= -1 &\text{ for }i=2, n_{arm}+2,2n_{arm}+2,...,(f_{arm}-1)N+2 \\
&M_{i,j} &= 0 &\text{ all others}
\end{align}
Examples
--------
Here we compare the FJC with FRC (costheta=0.1) of a 5 arm star in water.
We see the tranlational diffusion component at longer times (extrapolated to short times) and the faster relaxation
of internal modes at short times approaching the cumulant at shortest times.
::
import jscatter as js
import numpy as np
q= np.r_[0.01,0.1:2:0.2]
t = np.r_[0:1:0.02,1:20:1,20:400:5]
p = js.grace()
p.xaxis(label='t / ns')
sqt = js.dynamic.multiArmStarORZ(t, q, 5, 50, l=0.5, costheta = 0, reducedfriction=0.)
tt = t<20 # for cumulant
for c,sq in enumerate(sqt,1):
p.plot(sq.X,sq._Sqt,li=[1,1,c],sy=0,le=f'{sq.q:.2f} nm\S-1')
p.plot(sq.X, sq._Sqt_inf, li=[2, 1, c], sy=0, )
p.plot(sq.X[tt], sq._Sqt_cum[tt], li=[2, 1, 1], sy=0, )
# small
sqt = js.dynamic.multiArmStarORZ(t, q, 5, 50, l=0.5, costheta = 0.1, reducedfriction=0)
for c,sq in enumerate(sqt,1):
p.plot(sq.X,sq._Sqt,li=[3,1,c],sy=0)
p.yaxis(label='S(Q,t),S(Q,0)',scale='log',min=0.01,max=1)
p.legend(charsize=0.7)
p.title('5 arm star ORZ model: no HI and costheta=0, 0.1')
p.text('trans diffusion',x=30,y=0.16,rot=330)
p.text('costheta=0',x=300,y=0.031,rot=330)
p.text('costheta=0.1',x=300,y=0.022,rot=330)
p.text('cumulant diffusion',x=-20,y=0.3,rot=90)
# p.save(js.examples.imagepath+'/ORZ_Star.png',size=(1.5,1.5),dpi=300)
.. image:: ../../examples/images/ORZ_Star.png
:align: center
:width: 60 %
:alt: ORZ 5 arm star
We compare a flexible 10 arm star with a star that has bonds close to the core stretched.
We use a simple linear profile while in [1] a two step profile is used.
We observe that the stiff core increases tranlational diffusion (Dcm: 2.1 -> 1.58 nm²/ns) as the star gets larger
(Rg: 2.43 -> 3.99 nm). Additional the internal contribution increase in amplitude which might be a result
of the increased size as the arms are more extended.
::
import jscatter as js
import numpy as np
q= np.r_[0.01,0.1:2:0.2]
t = np.r_[0:1:0.02,1:20:1,20:100:5]
p = js.grace()
p.xaxis(label='t / ns')
sqt0 = js.dynamic.multiArmStarORZ(t, q, 10, 50, l=0.5, costheta = 0, reducedfriction=0.2)
for c,sq in enumerate(sqt0,1):
p.plot(sq.X,sq._Sqt,li=[1,1,c],sy=0,le=f'{sq.q:.2f} nm\\S-1')
# costheta is linear increasing from stretched center to free ends
sqt1 = js.dynamic.multiArmStarORZ(t, q, 5, 50, l=0.5, costheta = np.r_[0.7:0.1:49j], reducedfriction=0.2)
for c,sq in enumerate(sqt1,1):
p.plot(sq.X,sq._Sqt,li=[3,1,c],sy=0)
p.yaxis(label='S(Q,t),S(Q,0)',scale='log',min=0.01,max=1)
p.legend(charsize=0.7)
p.title('10 arm star ORZ model:')
p.subtitle('solid: costheta=0; broken costheta linear increasing')
p.text('costheta=0',x=300,y=0.031,rot=330)
p.text('costheta=0.1',x=300,y=0.022,rot=330)
# p.save(js.examples.imagepath+'/ORZ_10armStarblob.png',size=(1.5,1.5),dpi=300)
.. image:: ../../examples/images/ORZ_10armStarblob.png
:align: center
:width: 60 %
:alt: ORZ arm star blob model
References
----------
.. [1] Static and Dynamic Structure Factors for Star Polymers in θ Conditions.
Guenza, M. & Perico, A.
Macromolecules 26, 4196–4202 (1993). https://doi.org/10.1021/ma00068a020
.. [2] A Local Approach to the Dynamics of Star Polymers.
Guenza, M., Mormino, M. & Perico, A.
Macromolecules 24, 6168–6174 (1991). https://doi.org/10.1021/ma00023a018
.. [3] Theoretical basis for the Rouse‐Zimm model in polymer solution dynamics.
Zwanzig, R. The Journal of Chemical Physics 60, 2717–2720 (1974) https://doi.org/10.1063/1.1681433
.. [4] A hierarchy of models for the dynamics of polymer chains in dilute solution.
Perico, A., Ganazzoli, F. & Allegra, G.
The Journal of Chemical Physics 87, 3677–3686 (1987). https://doi.org/10.1063/1.452966
.. [5] Rubinstein, M. & Colby, R. H. Polymer Physics. (OUP Oxford, 2003).
.. [6] Intramolecular relaxation dynamics in semiflexible dendrimers.
Kumar, A. & Biswas, P. Journal of Chemical Physics 134, (2011). https://doi.org/10.1063/1.3598336
"""
assert np.all(np.array(costheta)<1) & np.all(0<= np.array(costheta)), 'costheta should be 0 <= cos(theta) < 1 !'
N = f_arm * n_arm + 1
if fa is None:
fa = np.ones(N)
assert len(fa) == N, 'fa should be of length f_arm * n_arm + 1.'
fa = np.array(fa)
q = np.atleast_1d(q)
t0 = np.r_[0,t]
# bead friction coefficient
friction = 6 * np.pi * (viscosity*1e-3) * (l * 1e-9) # l in nm; viscosity in Pa*s
# bond rate constant brc, sometimes named W
brc = (3 * kb * T / (l*1e-9)**2 / friction) * 1e-9 # l in nm # brc is 1/s => *1e-9 in 1/ns
reducedfriction = max(min(reducedfriction, 1), 0)
# create correction for diffusion
if Dcmfkt is not None:
if formel._getFuncCode(Dcmfkt):
# is already an interpolation function
Dcmfunktion = Dcmfkt
elif np.shape(Dcmfkt)[0] == 2:
Dcmfunktion = lambda qq: dA(Dcmfkt).interp(qq)
else:
raise TypeError('Shape of Dcmfkt is not 2xN!')
else:
# by default no correction
Dcmfunktion = lambda qq: 1.
# create bond matrix and transfer matrix for a star
A = _starStructuralMatrix(f_arm, n_arm, costheta)
# compute Eigenvectors and Eigenvalues for ORZ in reduced units
evals, evec, mu, loverR, Rg2_red, Rij2_red2, _ = solveOptimizedRouseZimm(A, reducedfriction)
if np.any(evals[1:] < 0):
raise UserWarning(f'There are {np.sum(evals[1:] < 0)} negative eigenvalues in ORZ solution. '
f'reducedfriction should be smaller.')
# >95% of computing time in this call for N=100 that already uses omp
# array dim [Q values x time values+2]
# first is sqt(0), last element is sum in equation 13 for initial cumulant, second last t=inf
# first eval is COM diffusion and large eval and large mu dont contribute much
lowev = np.r_[[False], (brc * evals[1:] < 15) | (1/mu[1:] > max(0.001/mu[1:]))]
sqt = fscatter.dynamic.sqtnonlinearpolymer(evec[:,1:], evals[1:], mu[1:], fa, q, t0, l, brc, loverR, tintern)
# [1] equ 24 for center of mass diffusion
Dcm_ORZ = kb * T / N / friction * 1e9 # in nm**2/ns
if reducedfriction>0:
Dcm_ORZ *= (1 + reducedfriction/N * np.sum(loverR))
results = dL()
for i, sqti in enumerate(sqt):
# calc initial slope diffusion as first cumulant/q**2 [1] equ 13
D_cum_ORZ = brc * l**2/3/sqti[0] * (np.sum(fa**2) + reducedfriction * sqti[-1])
if Dcm is None:
Dcm = Dcm_ORZ
results.append(np.c_[t0[1:],
np.exp(-q[i]**2 * Dcm * Dcmfunktion(q[i]) * t0[1:]) * sqti[1:-2]/sqti[0],
np.exp(-q[i]**2 * Dcm * Dcmfunktion(q[i]) * t0[1:]) * sqti[-2]/sqti[0],
np.exp(-q[i]**2 * (D_cum_ORZ - Dcm + Dcm * Dcmfunktion(q[i])) * t0[1:]) ].T)
results[-1].setColumnIndex(iey=None)
results[-1].columnname = 't;Sqt;Sqt_inf;Sqt_cum'
results[-1].viscosity = viscosity
results[-1].q = q[i]
results[-1].D_ORZcum = D_cum_ORZ
results[-1].D_ORZ = Dcm_ORZ
results[-1].Dcm = Dcm
results[-1].Fq = sqti[0] / np.sum(fa)**2 # form factor t=0
results[-1].Fq_inf = sqti[-2] / np.sum(fa)**2 # form factor t=inf
results[-1].beadfriction = friction
results[-1].bondrateconstant = brc # in 1/ns
results[-1].moderelaxationtimes = 1 / (brc * evals[lowev]) + tintern # [4] equ 2.17
results[-1].mode_rmsd = l / mu[lowev]**0.5 # [4] equ 2.19
results[-1].reducedfriction = reducedfriction
results[-1].costheta = costheta
results[-1].eigenvalues = evals
results[-1].mu = mu
results[-1].l = l
results[-1].tintern = tintern
results[-1].N = evec.shape[0]
results[-1].Rg = l * Rg2_red**0.5
results[-1].Rg_red = Rg2_red**0.5
results[-1].modelname = inspect.currentframe().f_code.co_name
if len(results) == 1:
return results[0]
return results
@formel.memoize(maxsize=4)
def _ringStructuralMatrix(N):
"""
Create for a ring the structural matrix for uncorrelated bonds.
This part is in reduced units that bond length is not used.
"""
# see references in main function
A = np.diag([2.]*N) + np.diag([-1.]*(N-1),-1) + np.diag([-1.]*(N-1),1)
A[0,-1] = -1
A[-1,0 ] = -1
return A
[docs]
def ringChainORZ(t, q, N=100, l=1, fa=None, Dcm=None, Dcmfkt=None,
viscosity=1, T=293, reducedfriction=0.25, tintern=0):
r"""
Dynamics of a ring polymer using optimized Rouse-Zimm approximation (ORZ).
The ring dynamic structure factor is calculated in analogy to the star described by Guenza [1]_.
The ring chain dynamics is decribed within the more general optimized Rouse-Zimm (ORZ) approximation
introduced by Bixon and Zwanzig [3]_ [4]_ . See :py:func:`solveOptimizedRouseZimm`.
We extend this here using bead scattering lengths for partial matching and allow individual bond correlations.
To speedup fitting use the memoize function as described in :py:func:`~.formel.memoize`.
Parameters
----------
t : array
timepoints in units ns.
q : array
Scattering vectors in units 1/nm.
N : int
Number of beads .
l : float, default = 0.38 nm (amino acid)
Bond length or Kuhn length in units nm.
fa : None, list of float
Scattering length of bead/monomer :math:`fa_i`. Can be used to match parts of the star to the solvent.
- None : Equal 1 for all beads.
- list: length N as scattering length of beads in sequence
Dcm : float
Center of mass diffusion in units nm**2/ns.
If `None` the calculated value D_ORZ is used.
Dcmfkt : array 2xN, function
Function f(q) or array with [qi, f(qi) ] as correction for Dcm like Diff = Dcm*f(q).
e.g. for inclusion of structure factor or hydrodynamic function with f(q)=H(Q)/S(q).
Missing values are interpolated.
viscosity : float, default=1 (H2O@20C )
:math:`\eta` in units cPoise=mPa*s e.g. water :math:`visc(T=293 K) =1 mPas`
T : float, default 273+20
Temperature in Kelvin.
reducedfriction : float, default 0.25
Reduced friction :math:`\zeta_r = \zeta/6\pi\eta_sl` in hydrodynamic tensor H.
- =0 : free draining limit, no HI, Rouse dynamics.
- >0 : with HI typically <0.5. For h>0.5 we get negative eigenvalues.
See :py:func:`solveOptimizedRouseZimm`.
During fits use ``limits(reducedfriction=[0,0.43,0.,0.5])`` to avoid singular matrices.
tintern : float>0, default 0
Relaxation time due to internal friction between neighboring beads in ns.
Returns
-------
sqt : dataList
Intermediate scattering function of a star for given q.
- [times; Sqt; Sqt only diffusion; Sqt cumulant]
- columnname = 't;Sqt;Sqt_inf;Sqt_cum'
- .q : scattering vector
- .D_ORZcum : diffusion coefficinet in initial slope (cumulant)
- .D_ORZ : translational diffusion ORZ model
- .Dcm : used center of friction/mass diffusion
- .Fq : form factor
- .Fq_inf : form factor t=inf
- .beadfriction : bead friction
- .bondrateconstant : bondrateconstant in 1/ns
- .moderelaxationtimes : used moderelaxationtimes in ns [4] equ 2.17
- .mode_rmsd : used mode rmsd as math:`l/\mu^{0.5} # [4] equ 2.19
- .reducedfriction : reducedfriction
- .eigenvalues : all evals
- .mu : all mu
- .l : bondlength l
- .N : Number of beads
- .Rg : radius of gyration in nm
- .Rg_red : reduced radius of gyration
Notes
-----
See :py:func:`solveOptimizedRouseZimm` for a description of the ORZ model with respective parameters and
the dynamic structure factor :math:`S(q,t)/S(q,0)` .
Here we use a ring chain with N beads of uncorelated beads with `costheta=0` .
The structural matrix has diagonal elements, :math:`A_{ii}=2` and :math:`A_{i\neq j}=-1` if the ith and jth monomers
are connected to each other or zero otherwise.
Examples
--------
Here we examine how HI changes dynamics.
::
import jscatter as js
import numpy as np
q= np.r_[0.01,0.1:2:0.2]
t = np.r_[0:1:0.02,1:20:1,20:100:5]
p = js.grace()
p.xaxis(label='t / ns')
sqt = js.dynamic.ringChainORZ(t, q, 100, l=0.5, reducedfriction=0)
for c,sq in enumerate(sqt,1):
p.plot(sq.X,sq._Sqt,li=[1,1,c],sy=0,le=f'{sq.q:.2f} nm\\S-1')
sqt = js.dynamic.ringChainORZ(t, q, 100, l=0.5, reducedfriction=0.05)
for c,sq in enumerate(sqt,1):
p.plot(sq.X,sq._Sqt,li=[3,1,c],sy=0,le='')
p.yaxis(label='S(Q,t),S(Q,0)',scale='log',min=0.1,max=1)
p.legend(charsize=0.7)
p.title('rings with and without HI ')
p.subtitle('solid: no HI; broken: with HI')
# p.save(js.examples.imagepath+'/ORZ_ringHI.png',size=(1.5,1.5),dpi=300)
.. image:: ../../examples/images/ORZ_ringHI.png
:align: center
:width: 60 %
:alt: ORZ eigenvalue and more
References
----------
.. [1] Static and Dynamic Structure Factors for Star Polymers in θ Conditions.
Guenza, M. & Perico, A.
Macromolecules 26, 4196–4202 (1993). https://doi.org/10.1021/ma00068a020
.. [2] A Local Approach to the Dynamics of Star Polymers.
Guenza, M., Mormino, M. & Perico, A.
Macromolecules 24, 6168–6174 (1991). https://doi.org/10.1021/ma00023a018
.. [3] Theoretical basis for the Rouse‐Zimm model in polymer solution dynamics.
Zwanzig, R. The Journal of Chemical Physics 60, 2717–2720 (1974) https://doi.org/10.1063/1.1681433
.. [4] A hierarchy of models for the dynamics of polymer chains in dilute solution.
Perico, A., Ganazzoli, F. & Allegra, G.
The Journal of Chemical Physics 87, 3677–3686 (1987). https://doi.org/10.1063/1.452966
.. [5] Rubinstein, M. & Colby, R. H. Polymer Physics. (OUP Oxford, 2003).
.. [6] Intramolecular relaxation dynamics in semiflexible dendrimers.
Kumar, A. & Biswas, P. Journal of Chemical Physics 134, (2011). https://doi.org/10.1063/1.3598336
"""
if fa is None:
fa = np.ones(N)
assert len(fa) == N, 'fa should be of length N.'
fa = np.array(fa)
q = np.atleast_1d(q)
t0 = np.r_[0,t]
# bead friction coefficient
friction = 6 * np.pi * (viscosity*1e-3) * (l * 1e-9) # l in nm; viscosity in Pa*s
# bond rate constant brc, sometimes named W
brc = (3 * kb * T / (l*1e-9)**2 / friction) * 1e-9 # l in nm # brc is 1/s => *1e-9 in 1/ns
reducedfriction = max(min(reducedfriction, 1), 0)
# create correction for diffusion
if Dcmfkt is not None:
if formel._getFuncCode(Dcmfkt):
# is already an interpolation function
Dcmfunktion = Dcmfkt
elif np.shape(Dcmfkt)[0] == 2:
Dcmfunktion = lambda qq: dA(Dcmfkt).interp(qq)
else:
raise TypeError('Shape of Dcmfkt is not 2xN!')
else:
# by default no correction
Dcmfunktion = lambda qq: 1.
# create bond matrix and transfer matrix for a star
A = _ringStructuralMatrix(N)
# compute Eigenvectors and Eigenvalues for ORZ in reduced units
evals, evec, mu, loverR, Rg2_red, Rij2_red2, _ = solveOptimizedRouseZimm(A, reducedfriction)
if np.any(evals[1:] < 0):
raise UserWarning(f'There are {np.sum(evals[1:]<0)} negative eigenvalues in ORZ solution. '
f'reducedfriction should be smaller.')
# >95% of computing time in this call for N=100 that already uses omp
# array dim [Q values x time values+2]
# first is sqt(0), last element is sum in equation 13 for initial cumulant, second last t=inf
# first eval is COM diffusion and large eval and large mu dont contribute much
lowev = np.r_[[False], (brc * evals[1:] < 15) | (1/mu[1:] > max(0.001/mu[1:]))]
sqt = fscatter.dynamic.sqtnonlinearpolymer(evec[:,1:], evals[1:], mu[1:], fa, q, t0, l, brc, loverR, tintern)
# [1] equ 24 for center of mass diffusion
Dcm_ORZ = kb * T / N / friction * 1e9 # in nm**2/ns
if reducedfriction>0:
Dcm_ORZ *= (1 + reducedfriction/N * np.sum(loverR))
results = dL()
for i, sqti in enumerate(sqt):
# calc initial slope diffusion as first cumulant/q**2 [1] equ 13
D_cum_ORZ = brc * l**2/3/sqti[0] * (np.sum(fa**2) + reducedfriction * sqti[-1])
if Dcm is None:
Dcm = Dcm_ORZ
results.append(np.c_[t0[1:],
np.exp(-q[i]**2 * Dcm * Dcmfunktion(q[i]) * t0[1:]) * sqti[1:-2]/sqti[0],
np.exp(-q[i]**2 * Dcm * Dcmfunktion(q[i]) * t0[1:]) * sqti[-2]/sqti[0],
np.exp(-q[i]**2 * (D_cum_ORZ - Dcm + Dcm * Dcmfunktion(q[i])) * t0[1:]) ].T)
results[-1].setColumnIndex(iey=None)
results[-1].columnname = 't;Sqt;Sqt_inf;Sqt_cum'
results[-1].viscosity = viscosity
results[-1].q = q[i]
results[-1].D_ORZcum = D_cum_ORZ
results[-1].D_ORZ = Dcm_ORZ
results[-1].Dcm = Dcm
results[-1].Fq = sqti[0] / np.sum(fa)**2 # form factor t=0
results[-1].Fq_inf = sqti[-2] / np.sum(fa)**2 # form factor t=inf
results[-1].beadfriction = friction
results[-1].bondrateconstant = brc # in 1/ns
results[-1].moderelaxationtimes = 1 / (brc * evals[lowev]) + tintern # [4] equ 2.17
results[-1].mode_rmsd = l / mu[lowev]**0.5 # [4] equ 2.19
results[-1].reducedfriction = reducedfriction
results[-1].eigenvalues = evals
results[-1].mu = mu
results[-1].l = l
results[-1].tintern = tintern
results[-1].N = evec.shape[0]
results[-1].Rg = l * Rg2_red**0.5
results[-1].Rg_red = Rg2_red**0.5
results[-1].modelname = inspect.currentframe().f_code.co_name
if len(results) == 1:
return results[0]
return results