Loading pysen/physics/nscatt.py +172 −58 Original line number Diff line number Diff line """ ================================== NSE scattering evaluator for a Rouse polymer chain ================================== ==================================================== ----------------------------------------------------------------------- Original autho: M. Monkenbusch Juelich Center of Neutron Science Forschungszentrum Juelich Original author: M. Monkenbusch, Juelich Center of Neutron Science Forschungszentrum Juelich m.monkenbusch@fz-juelich.de Fortran to Python translation: Piotr Zolnierczuk Piotr Zolnierczuk, Oak Ridge National Laboratory, zolnierczukp@ornl.gov ----------------------------------------------------------------------- ----------------------------------------------------------- This is a little toolbox to estimate the effect of contrast and H/D content on the dynamics data of a polymer melt in the Rouse regime as obtained by NSE spectroscopy. Loading @@ -26,13 +22,11 @@ given e.g. in the book of Doi-Edwards. But it will give a first impression on scattering contributions etc. In particular: consideration concerning the roles of coherent an incoherent scattering in polymer samples In particular: consideration concerning the roles of coherent and incoherent scattering in polymer samples please observe: ----------------------------------------------------------- units for units Q --> 1/Angstroem lambda --> Angstroem Loading @@ -57,16 +51,20 @@ t --> nanoseconds """ import numpy as np from numpy import pi, exp, cos, expm1, sqrt, mgrid, where #pylint: disable=no-name-in-module from numpy import pi, cos, exp, sqrt, mgrid, atleast_1d, asarray, sum as npsum, where, expm1 from scipy.special import erfc from scipy.integrate import quad from scipy.constants import Boltzmann as k_B from scipy.constants import (Boltzmann as k_B, angstrom as ANGSTROM) #we use a bunch of physics notations #pylint: disable=invalid-name _debye_cut = 1e-5 # cut-off value for f_debye _rouse_cut = 1e-33 # cut-off value for skernel_rouse _rouse_eps = 1e-6 # relative eps for Rouse kernel integration def f_debye(x): """Debye function f(y) where y=x^2=(kR)^2""" Loading @@ -75,32 +73,23 @@ def f_debye(x): return 2/y**2*(expm1(-y)+y) #pylint: disable=invalid-unary-operand-type def rouse_parameters(Wl4, Re, N, T=300.0): """ Calculate Rouse parameters using Wl4, Re and N """ l = Re/sqrt(N) # segment length kT = k_B * T * 100.0 # kg*m^2/s^2 -> kg * A^2/ns^2 zeta = 3*kT*l**2/Wl4 # friction factor W = 3*kT/(zeta*l**2) # Rouse rate in 1/ns Dr = kT/(N*zeta) # diffusion constant tR = (N/pi)**2/W # Rouse time Rg = Re/sqrt(6) # radius of gyration return dict(W=W, l=l, zeta_0=zeta, Rg=Rg, D=Dr, tau=tR, Wl4=Wl4, Re=Re, N=N, T=T) def rouse_sqt(t, q, W, Re, Dr=0, N=100, **kwargs): def _rouse_sqt(t, q, W, Re, **kwargs): """ Rouse expression for a chain of finite length: args:: t - time in ns q - momentum transfer in A^(-1) W - Rouse rate in 1/ns Re - end-to-end distance of the polymer molecule Re - end-to-end distance of the polymer molecule [A] kwargs:: Dr - center of mass diffusion constant in A**2/ns N - number of chain segments pmin - minimum p (default 1) pmax - maximum p (default N) ninf - if True, use different efac (see below) returns: Sq - S(Q) Sqt - S(Q,t) Loading @@ -108,29 +97,55 @@ def rouse_sqt(t, q, W, Re, Dr=0, N=100, **kwargs): Based on a Fortran subroutine nrousep by M. Monkenbusch (JCNS) Fortran -> Python translation by P. Zolnierczuk (JCNS) """ #we use a bunch of physics notations #pylint: disable=too-many-locals Dr = kwargs.pop('Dr', 0) N = kwargs.pop('N', 100) pmin = kwargs.pop('pmin', None) pmax = kwargs.pop('pmax', None) ninf = kwargs.pop('ninf', False) if N <= 0: raise RuntimeError('Number of chain segments must be positive') pmin = max(kwargs.get('pmin', 1), 1) pmax = min(kwargs.get('pmax', N), N) + 1 pmin = max(1,pmin or 1) pmax = min(N,pmax or N) + 1 l = Re/sqrt(N) # calculate segment length fac = 2/3*(Re*q/pi)**2 # # ---- Do the sums ----- PP, NN, MM = mgrid[pmin:pmax, 0:N, 0:N] # exponent for S(Q) and S(Q,t) anm = exp(-(q**2*abs(NN[0]-MM[0])/6*l**2)) # exponent for S(Q,t) cosf = exp(-fac*np.sum(cos(pi*(NN+1)*PP/N)*cos(pi*(MM+1)*PP/N)/PP**2 *(1-exp(-2*W*(1-cos(pi*PP/N))*t)), axis=(0,))) cospp = cos(pi*(NN+1)*PP/N)*cos(pi*(MM+1)*PP/N)/PP**2 if ninf: # this is for "infinite" N efac = (PP*pi/N)**2 #tR = (N/pi)**2/W else: # this is for "finite" N efac = 2*(1-cos(pi*PP/N)) cosf = exp(-fac*npsum(cospp*(1-exp(-efac*W*t)), axis=(0,))) Sq = np.sum(anm) Sqt = exp(-q**2*Dr*t) * np.sum(anm*cosf) Sq = npsum(anm) Sqt = exp(-q**2*Dr*t) * npsum(anm*cosf) return Sq/N, Sqt/N def rouse_sqt(t, q, W, Re, **kwargs): "vectorize in tau" res = [ _rouse_sqt(_t, q, W, Re, **kwargs) for _t in atleast_1d(t) ] return asarray(res) rouse_sqt.__doc__ = _rouse_sqt.__doc__ # ################################################################################### # WARNING THIS IS OBSOLETE APPROXIMATION FOR QR>>1 # THIS IS APPROXIMATION FOR QR>>1 # one should use rouse_sqt # Rouse kernel $h(u) = \frac{2}{\pi} \int_0^\infty \cos(ux)/x^2 (1-e^{-x^2}) dx$ # $h(u) = \frac{2}{\sqrt{pi}} \exp{-(u/2)^2} - u erfc(u/2) $ def h_rouse_hiq(u): Loading @@ -147,10 +162,109 @@ def skernel_rouse_hiq(u, xot): return exp(-(u+xot*h_rouse_hiq(u/xot))) # $ S(Q,t)/S(Q,0) = \int_0^\infty dx exp(-u + \sqrt(omegat) def sqt_rouse_hiq(xomegat): def rouse_sqt_hiq(omegat): "Rouse dynamic structure factor" if np.isscalar(xomegat): return quad(skernel_rouse_hiq, 0, np.inf, epsrel=_rouse_eps,args=(xomegat,)) res = [ quad(skernel_rouse_hiq, 0, np.inf, epsrel=_rouse_eps,args=(_xot, )) for _xot in xomegat ] return np.asarray(res) res = [ quad(skernel_rouse_hiq, 0, np.inf, epsrel=_rouse_eps, args=(_x,)) for _x in atleast_1d(omegat) ] return asarray(res) class RouseModel: """Rouse parameters units: - time ns - length: A (RE, Rg) - temperature: K - zeta: N ns/A """ def __init__(self, W, l, N, T=300): self.W = W # elemetary Rouse rate [1/ns] self.l = l # segment length [A] self.N = N # number of segments self.kT = k_B*T # kT def __str__(self): return ", ".join([f"{k}={v:g}" for k,v in self.to_dict().items()]) def units(self): "units" return "units=A,ns,K, zeta=N*ns/A" def to_dict(self): "convert to dictionary" return {'W' :self.W , 'l' :self.l, 'N' :self.N, 'kT':self.kT, 'RE' :self.RE, 'Rg':self.Rg, 'DR' :self.DR, 'tR':self.tR, 'Wl4' :self.Wl4, 'zeta': self.zeta} @classmethod def from_wl4(cls, Wl4, RE, N, T): "Calculate Rouse parameters using Wl4, Re and N" l = RE/sqrt(N) # segment length return RouseModel(Wl4/l**4,l,N,T) @property def RE(self): "end-to-end radius in A" return sqrt(self.N*self.l**2) @property def Rg(self): "Radius of gyration in A" return self.RE/sqrt(6) @property def zeta(self): "zeta in N*ns/A" return 3*self.kT/(self.W*self.l**2)/ANGSTROM @property def DR(self): "Rouse Diffusion coefficient [A^2/ns]" return self.Wl4/(3*self.RE**2) @property def tR(self): "Rouse time [ns]" return (self.N/pi)**2/self.W @property def Wl4(self): "Wl^4 nm^4/ns" return self.W*self.l**4 def omegaR(self, q): "Rouse <variable>" return (q*self.l)**4*self.W/36 def sqt(self, t, q, **kwargs): """ Rouse expression for a chain of finite length: args:: t - time in ns q - momentum transfer in A^(-1) kwargs:: Dr - center of mass diffusion constant in A**2/ns pmin - minimum p (default 1) pmax - maximum p (default N) ninf - if True, use different efac (see below) returns: Sq - S(Q) Sqt - S(Q,t) """ return rouse_sqt(t, q, W=self.W, Re=self.RE, N=self.N, **kwargs) def sqt_hiq(self, t, q): """ Rouse expression for high Q. args:: t - time in ns q - momentum transfer in A^(-1) returns: Sqt - S(Q,t) """ u = sqrt(self.omegaR(q)*atleast_1d(t)) res = np.asarray(rouse_sqt_hiq(u)) return res # pysen/revision.py +2 −2 Original line number Diff line number Diff line Loading @@ -3,8 +3,8 @@ PySEN revision module """ import sys __version__ = "2.0" __release__ = "a8" __date__ = "Nov 14, 2024" __release__ = "a9" __date__ = "Dec 19, 2024" def version(full=False): "get pysen version number" Loading test/test_nscatt.py +104 −47 Original line number Diff line number Diff line Loading @@ -6,76 +6,78 @@ import unittest import numpy as np import scipy.integrate as integ import pysen.physics.nscatt as nscatt from numpy import sqrt, log10, pi from scipy.constants import Boltzmann as kB from pysen.physics.nscatt import (RouseModel, f_debye, h_rouse_hiq, rouse_sqt_hiq, rouse_sqt) class NeutronScatteringTestCase(unittest.TestCase): "Neutron Scattering Test Cases" def setUp(self): self.Wl4 = 39500.0 self.Re = 40.0 self.T = 300.0 self.N = 100 pass # values from Richter, Monkenbusch, Arbe, Colmenero (2005) # for PEP at 423K self.T = 423 # temperature [K] self.N = 100 # number of segments self.Wl4 = 0.95e4 # Rouse parameter Wl^4 [A^4/ns] self.RE = sqrt(self.N*56.5) # end-to-end radius [A] def tearDown(self): pass def test_f_debye(self): "test for Debye function" inp = np.array([ # x f(x) ( 0.0, 1.0000), ( 1.0, 0.7358), ( 2.0, 0.3773), (50.0, 0.0008), ]) fx = nscatt.f_debye(inp[:,0]) fx = f_debye(inp[:,0]) self.assertTrue(np.allclose(fx, inp[:,1], atol=1e-03)) # I0 = integ.quad(nscatt.f_debye, 0, np.infty) self.assertAlmostEqual(4/3*np.sqrt(np.pi), I0[0], 5, "f_debye integral (value)") I0 = integ.quad(f_debye, 0, np.infty) self.assertAlmostEqual(4/3*sqrt(pi), I0[0], 5, "f_debye integral (value)") self.assertAlmostEqual(0 , I0[1], 5, "f_debye integral (error)") def test_rouse_hkernel_hiq(self): def test_rouse_parameters(self): "test for Rouse model parameters" rpar = RouseModel.from_wl4(self.Wl4, self.RE, self.N, self.T) #print("\n",rpar) N = self.N l2 = self.RE**2/N # l^2 logW = log10(self.Wl4/(l2**2)) # W = Wl4/(l^2)^2 logDR = log10(self.Wl4/(3*self.RE**2)) # DR = Wl4/(3*RE^2) logtR = log10((self.RE**2/pi)**2/self.Wl4) # tR = (RE/pi)^2/Wl4 logZeta = log10(3*kB*self.T/self.Wl4*l2)+10 # +10 from m to A # self.assertAlmostEqual(rpar.kT/kB, self.T , 5, "RouseParameter T") self.assertAlmostEqual(rpar.N, self.N , 5, "RouseParameter N") self.assertAlmostEqual(rpar.Wl4, self.Wl4 , 5, "RouseParameter Wl4") self.assertAlmostEqual(rpar.RE, self.RE , 5, "RouseParameter RE") self.assertAlmostEqual(rpar.l, self.RE/sqrt(N) , 5, "RouseParameter ell") self.assertAlmostEqual(rpar.Rg, self.RE/sqrt(6) , 5, "RouseParameter Rg") self.assertAlmostEqual(log10(rpar.W), logW , 5, "RouseParameter W ") self.assertAlmostEqual(log10(rpar.DR), logDR , 5, "RouseParameter DR") self.assertAlmostEqual(log10(rpar.tR), logtR , 5, "RouseParameter tR") self.assertAlmostEqual(log10(rpar.zeta), logZeta , 5, "RouseParameter zeta") def test_rouse_kernel_hiq(self): "test for Rouse hiQ kernel" inp = np.array([ # u h(u) ( 0.0, 2/np.sqrt(np.pi)), ( 1.0, 0.3993), ( 2.0, 0.1005), (10.0, 0.0000), ]) hu = nscatt.h_rouse_hiq(inp[:,0]) hu = h_rouse_hiq(inp[:,0]) self.assertTrue(np.allclose(hu, inp[:,1], atol=1e-03)) def test_sqt_rouse_hiq(self): inp = np.array([ # omegat S(omegat) ( 0.0, 1.0000), ( 1.0, 0.6088), ( 2.0, 0.2620), ( 5.0, 0.0135), (10.0, 0.0000), ]) sqt = nscatt.sqt_rouse_hiq(inp[:,0])[:,0] self.assertTrue(np.allclose(sqt, inp[:,1], atol=1e-03), sqt) def test_rouse_parameters(self): r = nscatt.rouse_parameters(Wl4=self.Wl4, Re=self.Re, T=self.T, N=self.N) # input self.assertEqual(self.N , r['N'] , "Number of segments") self.assertAlmostEqual(self.Wl4 , r['Wl4'] , 5, "Wl4: Rouse rate") self.assertAlmostEqual(self.Re , r['Re'] , 5, "end to end length") self.assertAlmostEqual(self.T , r['T'] , 5, "temperature") # results self.assertAlmostEqual( 4.0 , r['l'] , 5, "l: segment length") self.assertAlmostEqual( 16.3299316 , r['Rg'] , 5, "radius of gyration") self.assertAlmostEqual( 5.033251 , r['zeta_0']*1e22, 5, "zeta: friction constant") self.assertAlmostEqual(154.296875 , r['W'] , 5, "W: Rouse rate") self.assertAlmostEqual( 6.566638 , r['tau'] , 5, "tR: Rouse time") self.assertAlmostEqual( 8.229167 , r['D'] , 5, "D: Rouse diffusion rate") def test_rouse_sqt(self): def test_rouse_sqt1(self): "test for Rouse Sqt (vs datreat)" # results from datreat # tau S(Q) S(Q,t) d = np.asarray( [[ 0.00000, 48.835216, 48.835216], Loading @@ -84,12 +86,67 @@ class NeutronScatteringTestCase(unittest.TestCase): [30.00000, 48.835216, 42.692559], [40.00000, 48.835216, 42.692053], [50.00000, 48.835216, 42.692025]]) W = nscatt.rouse_parameters(Wl4=self.Wl4, Re=self.Re, T=self.T, N=self.N)['W'] q = 0.1 rmod = RouseModel.from_wl4(Wl4=39500.0, RE=40.0, N=100, T=300) # first loop direct rouse_sqt for tau, Sq, Sqt in d: act = rouse_sqt(tau, q , W=rmod.W, Re=rmod.RE, N=rmod.N, pmin=1, pmax=rmod.N) self.assertAlmostEqual(Sq, act[0,0], 5, "S(Q)") self.assertAlmostEqual(Sqt, act[0,1], 5, "S(Q,t)") # second loop direct model method for tau, Sq, Sqt in d: act_Sq, act_Sqt = nscatt.rouse_sqt(tau, q , W, self.Re, N=self.N, pmin=1, pmax=self.N) self.assertAlmostEqual(Sq, act_Sq, 4, "S(Q)") self.assertAlmostEqual(Sqt, act_Sqt, 5, "S(Q,t)") act = rmod.sqt(tau, q , pmin=1, pmax=rmod.N) self.assertAlmostEqual(Sq, act[0,0], 5, "S(Q)") self.assertAlmostEqual(Sqt, act[0,1], 5, "S(Q,t)") def test_rouse_sqt2(self): "test for Rouse Sqt" # results for PEP@423 # S(Q) S(Q,t) data = [[5.00501524, 5.00501524], [5.00501524, 2.16188872], [5.00501524, 1.34832783], [5.00501524, 0.93827628], [5.00501524, 0.69473812], [5.00501524, 0.53632341]] expected = np.asarray(data) rmod = RouseModel.from_wl4(self.Wl4, self.RE, self.N, self.T) actual = rmod.sqt([0,5,10,15,20,25], 2*pi/rmod.Rg , pmax=rmod.N) self.assertTrue(np.allclose(actual, expected, atol=1e-10), "rouse_sqt2") def test_sqt_rouse_hiq1(self): "test for Rouse Sqt (high Q)" inp = np.array([ # omegat S(omegat) ( 0.0, 1.0000), ( 1.0, 0.6088), ( 2.0, 0.2620), ( 5.0, 0.0135), (10.0, 0.0000), ]) sqt = rouse_sqt_hiq(inp[:,0])[:,0] self.assertTrue(np.allclose(sqt, inp[:,1], atol=1e-03), "sqt_hiq1") def test_sqt_rouse_hiq2(self): "test for Rouse Sqt (high Q)" inp = np.array([ # omegat S(omegat) ( 0.0, 1.0000), ( 1.0, 0.6088), ( 2.0, 0.2620), ( 5.0, 0.0135), (10.0, 0.0000), ]) u, expected = inp[:,0], inp[:,1] rmod = RouseModel.from_wl4(self.Wl4, self.RE, self.N, self.T) q=2*pi/rmod.Rg taus = u**2/rmod.omegaR(q) actual = rmod.sqt_hiq(taus, q)[:,0] self.assertTrue(np.allclose(actual, expected, atol=1e-4), "rouse_hiq2") if __name__ == "__main__": unittest.main(exit=False, verbosity=2) Loading
pysen/physics/nscatt.py +172 −58 Original line number Diff line number Diff line """ ================================== NSE scattering evaluator for a Rouse polymer chain ================================== ==================================================== ----------------------------------------------------------------------- Original autho: M. Monkenbusch Juelich Center of Neutron Science Forschungszentrum Juelich Original author: M. Monkenbusch, Juelich Center of Neutron Science Forschungszentrum Juelich m.monkenbusch@fz-juelich.de Fortran to Python translation: Piotr Zolnierczuk Piotr Zolnierczuk, Oak Ridge National Laboratory, zolnierczukp@ornl.gov ----------------------------------------------------------------------- ----------------------------------------------------------- This is a little toolbox to estimate the effect of contrast and H/D content on the dynamics data of a polymer melt in the Rouse regime as obtained by NSE spectroscopy. Loading @@ -26,13 +22,11 @@ given e.g. in the book of Doi-Edwards. But it will give a first impression on scattering contributions etc. In particular: consideration concerning the roles of coherent an incoherent scattering in polymer samples In particular: consideration concerning the roles of coherent and incoherent scattering in polymer samples please observe: ----------------------------------------------------------- units for units Q --> 1/Angstroem lambda --> Angstroem Loading @@ -57,16 +51,20 @@ t --> nanoseconds """ import numpy as np from numpy import pi, exp, cos, expm1, sqrt, mgrid, where #pylint: disable=no-name-in-module from numpy import pi, cos, exp, sqrt, mgrid, atleast_1d, asarray, sum as npsum, where, expm1 from scipy.special import erfc from scipy.integrate import quad from scipy.constants import Boltzmann as k_B from scipy.constants import (Boltzmann as k_B, angstrom as ANGSTROM) #we use a bunch of physics notations #pylint: disable=invalid-name _debye_cut = 1e-5 # cut-off value for f_debye _rouse_cut = 1e-33 # cut-off value for skernel_rouse _rouse_eps = 1e-6 # relative eps for Rouse kernel integration def f_debye(x): """Debye function f(y) where y=x^2=(kR)^2""" Loading @@ -75,32 +73,23 @@ def f_debye(x): return 2/y**2*(expm1(-y)+y) #pylint: disable=invalid-unary-operand-type def rouse_parameters(Wl4, Re, N, T=300.0): """ Calculate Rouse parameters using Wl4, Re and N """ l = Re/sqrt(N) # segment length kT = k_B * T * 100.0 # kg*m^2/s^2 -> kg * A^2/ns^2 zeta = 3*kT*l**2/Wl4 # friction factor W = 3*kT/(zeta*l**2) # Rouse rate in 1/ns Dr = kT/(N*zeta) # diffusion constant tR = (N/pi)**2/W # Rouse time Rg = Re/sqrt(6) # radius of gyration return dict(W=W, l=l, zeta_0=zeta, Rg=Rg, D=Dr, tau=tR, Wl4=Wl4, Re=Re, N=N, T=T) def rouse_sqt(t, q, W, Re, Dr=0, N=100, **kwargs): def _rouse_sqt(t, q, W, Re, **kwargs): """ Rouse expression for a chain of finite length: args:: t - time in ns q - momentum transfer in A^(-1) W - Rouse rate in 1/ns Re - end-to-end distance of the polymer molecule Re - end-to-end distance of the polymer molecule [A] kwargs:: Dr - center of mass diffusion constant in A**2/ns N - number of chain segments pmin - minimum p (default 1) pmax - maximum p (default N) ninf - if True, use different efac (see below) returns: Sq - S(Q) Sqt - S(Q,t) Loading @@ -108,29 +97,55 @@ def rouse_sqt(t, q, W, Re, Dr=0, N=100, **kwargs): Based on a Fortran subroutine nrousep by M. Monkenbusch (JCNS) Fortran -> Python translation by P. Zolnierczuk (JCNS) """ #we use a bunch of physics notations #pylint: disable=too-many-locals Dr = kwargs.pop('Dr', 0) N = kwargs.pop('N', 100) pmin = kwargs.pop('pmin', None) pmax = kwargs.pop('pmax', None) ninf = kwargs.pop('ninf', False) if N <= 0: raise RuntimeError('Number of chain segments must be positive') pmin = max(kwargs.get('pmin', 1), 1) pmax = min(kwargs.get('pmax', N), N) + 1 pmin = max(1,pmin or 1) pmax = min(N,pmax or N) + 1 l = Re/sqrt(N) # calculate segment length fac = 2/3*(Re*q/pi)**2 # # ---- Do the sums ----- PP, NN, MM = mgrid[pmin:pmax, 0:N, 0:N] # exponent for S(Q) and S(Q,t) anm = exp(-(q**2*abs(NN[0]-MM[0])/6*l**2)) # exponent for S(Q,t) cosf = exp(-fac*np.sum(cos(pi*(NN+1)*PP/N)*cos(pi*(MM+1)*PP/N)/PP**2 *(1-exp(-2*W*(1-cos(pi*PP/N))*t)), axis=(0,))) cospp = cos(pi*(NN+1)*PP/N)*cos(pi*(MM+1)*PP/N)/PP**2 if ninf: # this is for "infinite" N efac = (PP*pi/N)**2 #tR = (N/pi)**2/W else: # this is for "finite" N efac = 2*(1-cos(pi*PP/N)) cosf = exp(-fac*npsum(cospp*(1-exp(-efac*W*t)), axis=(0,))) Sq = np.sum(anm) Sqt = exp(-q**2*Dr*t) * np.sum(anm*cosf) Sq = npsum(anm) Sqt = exp(-q**2*Dr*t) * npsum(anm*cosf) return Sq/N, Sqt/N def rouse_sqt(t, q, W, Re, **kwargs): "vectorize in tau" res = [ _rouse_sqt(_t, q, W, Re, **kwargs) for _t in atleast_1d(t) ] return asarray(res) rouse_sqt.__doc__ = _rouse_sqt.__doc__ # ################################################################################### # WARNING THIS IS OBSOLETE APPROXIMATION FOR QR>>1 # THIS IS APPROXIMATION FOR QR>>1 # one should use rouse_sqt # Rouse kernel $h(u) = \frac{2}{\pi} \int_0^\infty \cos(ux)/x^2 (1-e^{-x^2}) dx$ # $h(u) = \frac{2}{\sqrt{pi}} \exp{-(u/2)^2} - u erfc(u/2) $ def h_rouse_hiq(u): Loading @@ -147,10 +162,109 @@ def skernel_rouse_hiq(u, xot): return exp(-(u+xot*h_rouse_hiq(u/xot))) # $ S(Q,t)/S(Q,0) = \int_0^\infty dx exp(-u + \sqrt(omegat) def sqt_rouse_hiq(xomegat): def rouse_sqt_hiq(omegat): "Rouse dynamic structure factor" if np.isscalar(xomegat): return quad(skernel_rouse_hiq, 0, np.inf, epsrel=_rouse_eps,args=(xomegat,)) res = [ quad(skernel_rouse_hiq, 0, np.inf, epsrel=_rouse_eps,args=(_xot, )) for _xot in xomegat ] return np.asarray(res) res = [ quad(skernel_rouse_hiq, 0, np.inf, epsrel=_rouse_eps, args=(_x,)) for _x in atleast_1d(omegat) ] return asarray(res) class RouseModel: """Rouse parameters units: - time ns - length: A (RE, Rg) - temperature: K - zeta: N ns/A """ def __init__(self, W, l, N, T=300): self.W = W # elemetary Rouse rate [1/ns] self.l = l # segment length [A] self.N = N # number of segments self.kT = k_B*T # kT def __str__(self): return ", ".join([f"{k}={v:g}" for k,v in self.to_dict().items()]) def units(self): "units" return "units=A,ns,K, zeta=N*ns/A" def to_dict(self): "convert to dictionary" return {'W' :self.W , 'l' :self.l, 'N' :self.N, 'kT':self.kT, 'RE' :self.RE, 'Rg':self.Rg, 'DR' :self.DR, 'tR':self.tR, 'Wl4' :self.Wl4, 'zeta': self.zeta} @classmethod def from_wl4(cls, Wl4, RE, N, T): "Calculate Rouse parameters using Wl4, Re and N" l = RE/sqrt(N) # segment length return RouseModel(Wl4/l**4,l,N,T) @property def RE(self): "end-to-end radius in A" return sqrt(self.N*self.l**2) @property def Rg(self): "Radius of gyration in A" return self.RE/sqrt(6) @property def zeta(self): "zeta in N*ns/A" return 3*self.kT/(self.W*self.l**2)/ANGSTROM @property def DR(self): "Rouse Diffusion coefficient [A^2/ns]" return self.Wl4/(3*self.RE**2) @property def tR(self): "Rouse time [ns]" return (self.N/pi)**2/self.W @property def Wl4(self): "Wl^4 nm^4/ns" return self.W*self.l**4 def omegaR(self, q): "Rouse <variable>" return (q*self.l)**4*self.W/36 def sqt(self, t, q, **kwargs): """ Rouse expression for a chain of finite length: args:: t - time in ns q - momentum transfer in A^(-1) kwargs:: Dr - center of mass diffusion constant in A**2/ns pmin - minimum p (default 1) pmax - maximum p (default N) ninf - if True, use different efac (see below) returns: Sq - S(Q) Sqt - S(Q,t) """ return rouse_sqt(t, q, W=self.W, Re=self.RE, N=self.N, **kwargs) def sqt_hiq(self, t, q): """ Rouse expression for high Q. args:: t - time in ns q - momentum transfer in A^(-1) returns: Sqt - S(Q,t) """ u = sqrt(self.omegaR(q)*atleast_1d(t)) res = np.asarray(rouse_sqt_hiq(u)) return res #
pysen/revision.py +2 −2 Original line number Diff line number Diff line Loading @@ -3,8 +3,8 @@ PySEN revision module """ import sys __version__ = "2.0" __release__ = "a8" __date__ = "Nov 14, 2024" __release__ = "a9" __date__ = "Dec 19, 2024" def version(full=False): "get pysen version number" Loading
test/test_nscatt.py +104 −47 Original line number Diff line number Diff line Loading @@ -6,76 +6,78 @@ import unittest import numpy as np import scipy.integrate as integ import pysen.physics.nscatt as nscatt from numpy import sqrt, log10, pi from scipy.constants import Boltzmann as kB from pysen.physics.nscatt import (RouseModel, f_debye, h_rouse_hiq, rouse_sqt_hiq, rouse_sqt) class NeutronScatteringTestCase(unittest.TestCase): "Neutron Scattering Test Cases" def setUp(self): self.Wl4 = 39500.0 self.Re = 40.0 self.T = 300.0 self.N = 100 pass # values from Richter, Monkenbusch, Arbe, Colmenero (2005) # for PEP at 423K self.T = 423 # temperature [K] self.N = 100 # number of segments self.Wl4 = 0.95e4 # Rouse parameter Wl^4 [A^4/ns] self.RE = sqrt(self.N*56.5) # end-to-end radius [A] def tearDown(self): pass def test_f_debye(self): "test for Debye function" inp = np.array([ # x f(x) ( 0.0, 1.0000), ( 1.0, 0.7358), ( 2.0, 0.3773), (50.0, 0.0008), ]) fx = nscatt.f_debye(inp[:,0]) fx = f_debye(inp[:,0]) self.assertTrue(np.allclose(fx, inp[:,1], atol=1e-03)) # I0 = integ.quad(nscatt.f_debye, 0, np.infty) self.assertAlmostEqual(4/3*np.sqrt(np.pi), I0[0], 5, "f_debye integral (value)") I0 = integ.quad(f_debye, 0, np.infty) self.assertAlmostEqual(4/3*sqrt(pi), I0[0], 5, "f_debye integral (value)") self.assertAlmostEqual(0 , I0[1], 5, "f_debye integral (error)") def test_rouse_hkernel_hiq(self): def test_rouse_parameters(self): "test for Rouse model parameters" rpar = RouseModel.from_wl4(self.Wl4, self.RE, self.N, self.T) #print("\n",rpar) N = self.N l2 = self.RE**2/N # l^2 logW = log10(self.Wl4/(l2**2)) # W = Wl4/(l^2)^2 logDR = log10(self.Wl4/(3*self.RE**2)) # DR = Wl4/(3*RE^2) logtR = log10((self.RE**2/pi)**2/self.Wl4) # tR = (RE/pi)^2/Wl4 logZeta = log10(3*kB*self.T/self.Wl4*l2)+10 # +10 from m to A # self.assertAlmostEqual(rpar.kT/kB, self.T , 5, "RouseParameter T") self.assertAlmostEqual(rpar.N, self.N , 5, "RouseParameter N") self.assertAlmostEqual(rpar.Wl4, self.Wl4 , 5, "RouseParameter Wl4") self.assertAlmostEqual(rpar.RE, self.RE , 5, "RouseParameter RE") self.assertAlmostEqual(rpar.l, self.RE/sqrt(N) , 5, "RouseParameter ell") self.assertAlmostEqual(rpar.Rg, self.RE/sqrt(6) , 5, "RouseParameter Rg") self.assertAlmostEqual(log10(rpar.W), logW , 5, "RouseParameter W ") self.assertAlmostEqual(log10(rpar.DR), logDR , 5, "RouseParameter DR") self.assertAlmostEqual(log10(rpar.tR), logtR , 5, "RouseParameter tR") self.assertAlmostEqual(log10(rpar.zeta), logZeta , 5, "RouseParameter zeta") def test_rouse_kernel_hiq(self): "test for Rouse hiQ kernel" inp = np.array([ # u h(u) ( 0.0, 2/np.sqrt(np.pi)), ( 1.0, 0.3993), ( 2.0, 0.1005), (10.0, 0.0000), ]) hu = nscatt.h_rouse_hiq(inp[:,0]) hu = h_rouse_hiq(inp[:,0]) self.assertTrue(np.allclose(hu, inp[:,1], atol=1e-03)) def test_sqt_rouse_hiq(self): inp = np.array([ # omegat S(omegat) ( 0.0, 1.0000), ( 1.0, 0.6088), ( 2.0, 0.2620), ( 5.0, 0.0135), (10.0, 0.0000), ]) sqt = nscatt.sqt_rouse_hiq(inp[:,0])[:,0] self.assertTrue(np.allclose(sqt, inp[:,1], atol=1e-03), sqt) def test_rouse_parameters(self): r = nscatt.rouse_parameters(Wl4=self.Wl4, Re=self.Re, T=self.T, N=self.N) # input self.assertEqual(self.N , r['N'] , "Number of segments") self.assertAlmostEqual(self.Wl4 , r['Wl4'] , 5, "Wl4: Rouse rate") self.assertAlmostEqual(self.Re , r['Re'] , 5, "end to end length") self.assertAlmostEqual(self.T , r['T'] , 5, "temperature") # results self.assertAlmostEqual( 4.0 , r['l'] , 5, "l: segment length") self.assertAlmostEqual( 16.3299316 , r['Rg'] , 5, "radius of gyration") self.assertAlmostEqual( 5.033251 , r['zeta_0']*1e22, 5, "zeta: friction constant") self.assertAlmostEqual(154.296875 , r['W'] , 5, "W: Rouse rate") self.assertAlmostEqual( 6.566638 , r['tau'] , 5, "tR: Rouse time") self.assertAlmostEqual( 8.229167 , r['D'] , 5, "D: Rouse diffusion rate") def test_rouse_sqt(self): def test_rouse_sqt1(self): "test for Rouse Sqt (vs datreat)" # results from datreat # tau S(Q) S(Q,t) d = np.asarray( [[ 0.00000, 48.835216, 48.835216], Loading @@ -84,12 +86,67 @@ class NeutronScatteringTestCase(unittest.TestCase): [30.00000, 48.835216, 42.692559], [40.00000, 48.835216, 42.692053], [50.00000, 48.835216, 42.692025]]) W = nscatt.rouse_parameters(Wl4=self.Wl4, Re=self.Re, T=self.T, N=self.N)['W'] q = 0.1 rmod = RouseModel.from_wl4(Wl4=39500.0, RE=40.0, N=100, T=300) # first loop direct rouse_sqt for tau, Sq, Sqt in d: act = rouse_sqt(tau, q , W=rmod.W, Re=rmod.RE, N=rmod.N, pmin=1, pmax=rmod.N) self.assertAlmostEqual(Sq, act[0,0], 5, "S(Q)") self.assertAlmostEqual(Sqt, act[0,1], 5, "S(Q,t)") # second loop direct model method for tau, Sq, Sqt in d: act_Sq, act_Sqt = nscatt.rouse_sqt(tau, q , W, self.Re, N=self.N, pmin=1, pmax=self.N) self.assertAlmostEqual(Sq, act_Sq, 4, "S(Q)") self.assertAlmostEqual(Sqt, act_Sqt, 5, "S(Q,t)") act = rmod.sqt(tau, q , pmin=1, pmax=rmod.N) self.assertAlmostEqual(Sq, act[0,0], 5, "S(Q)") self.assertAlmostEqual(Sqt, act[0,1], 5, "S(Q,t)") def test_rouse_sqt2(self): "test for Rouse Sqt" # results for PEP@423 # S(Q) S(Q,t) data = [[5.00501524, 5.00501524], [5.00501524, 2.16188872], [5.00501524, 1.34832783], [5.00501524, 0.93827628], [5.00501524, 0.69473812], [5.00501524, 0.53632341]] expected = np.asarray(data) rmod = RouseModel.from_wl4(self.Wl4, self.RE, self.N, self.T) actual = rmod.sqt([0,5,10,15,20,25], 2*pi/rmod.Rg , pmax=rmod.N) self.assertTrue(np.allclose(actual, expected, atol=1e-10), "rouse_sqt2") def test_sqt_rouse_hiq1(self): "test for Rouse Sqt (high Q)" inp = np.array([ # omegat S(omegat) ( 0.0, 1.0000), ( 1.0, 0.6088), ( 2.0, 0.2620), ( 5.0, 0.0135), (10.0, 0.0000), ]) sqt = rouse_sqt_hiq(inp[:,0])[:,0] self.assertTrue(np.allclose(sqt, inp[:,1], atol=1e-03), "sqt_hiq1") def test_sqt_rouse_hiq2(self): "test for Rouse Sqt (high Q)" inp = np.array([ # omegat S(omegat) ( 0.0, 1.0000), ( 1.0, 0.6088), ( 2.0, 0.2620), ( 5.0, 0.0135), (10.0, 0.0000), ]) u, expected = inp[:,0], inp[:,1] rmod = RouseModel.from_wl4(self.Wl4, self.RE, self.N, self.T) q=2*pi/rmod.Rg taus = u**2/rmod.omegaR(q) actual = rmod.sqt_hiq(taus, q)[:,0] self.assertTrue(np.allclose(actual, expected, atol=1e-4), "rouse_hiq2") if __name__ == "__main__": unittest.main(exit=False, verbosity=2)
