more code in RouseModel (testing and added debye S(Q)) (d6928425) · Commits · Zolnierczuk, Piotr / pysen

pysen/physics/nscatt.py

+47 −32

Original line number	Diff line number	Diff line
		@@ -60,7 +60,7 @@ 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
		_debye_cut = 1e-6 # 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

		@@ -73,7 +73,44 @@ def f_debye(x):
		return 2/y*2(expm1(-y)+y) #pylint: disable=invalid-unary-operand-type


		# ###################################################################################
		# THIS IS APPROXIMATION FOR QR>>1
		# one should use rouse_sqt (below)
		# 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):
		"Rouse kernel h(u)"
		u2 = u/2
		return 2/sqrt(pi)exp(-(u2)2) - uerfc(u2)

		# exp(-(u+sqrt(omegat)*h_rouse(u/sqrt(omegat))))
		def skernel_rouse_hiq(u, xot):
		"Rouse S(Q,t) inegrand"
		#xot = sqrt(omegat)
		if xot<_rouse_cut:
		return exp(-u) # h(oo) = 0
		return exp(-(u+xot*h_rouse_hiq(u/xot)))

		# $ S(Q,t)/S(Q,0) = \int_0^\infty dx exp(-u + \sqrt(omegat)
		def rouse_sqt_hiq(omegat, accuracy=False):
		"""Rouse dynamic structure factor
		The original expression contains 12/(Q*l)^2 prefactor which cancels
		when getting S(Q,t)/S(Q,0)
		"""
		res = [ quad(skernel_rouse_hiq, 0, np.inf, epsrel=_rouse_eps, args=(_x,))
		for _x in atleast_1d(omegat) ]
		res = asarray(res)
		if accuracy:
		# return Sqt and accuracy
		return res[:,0], res[:,1]
		# return Sqt only
		return res[:,0]

		# ###################################################################################
		# Rouse model S(Q,t)
		# mode summation
		# ###################################################################################
		def _rouse_sqt(t, q, W, Re, **kwargs):
		"""
		Rouse expression for a chain of finite length:
		@@ -143,31 +180,6 @@ def rouse_sqt(t, q, W, Re, **kwargs):
		return asarray(res)
		rouse_sqt.__doc__ = _rouse_sqt.__doc__

		# ###################################################################################
		# 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):
		"Rouse kernel h(u)"
		u2 = u/2
		return 2/sqrt(pi)exp(-(u2)2) - uerfc(u2)

		# exp(-(u+sqrt(omegat)*h_rouse(u/sqrt(omegat))))
		def skernel_rouse_hiq(u, xot):
		"Rouse S(Q,t) inegrand"
		#xot = sqrt(omegat)
		if xot<_rouse_cut:
		return exp(-u) # h(oo) = 0
		return exp(-(u+xot*h_rouse_hiq(u/xot)))

		# $ S(Q,t)/S(Q,0) = \int_0^\infty dx exp(-u + \sqrt(omegat)
		def rouse_sqt_hiq(omegat):
		"Rouse dynamic structure factor"
		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
		@@ -231,9 +243,13 @@ class RouseModel:
		return self.Wself.l*4

		def omegaR(self, q):
		"Rouse <variable>"
		"Rouse <variable> also known as Gamma_q"
		return (qself.l)4self.W/36

		def sq(self, q):
		"Debye structure factor"
		return self.Nf_debye((self.Rgq)**2)

		def sqt(self, t, q, **kwargs):
		"""
		Rouse expression for a chain of finite length:
		@@ -264,7 +280,6 @@ class RouseModel:
		Sqt - S(Q,t)
		"""
		u = sqrt(self.omegaR(q)*atleast_1d(t))
		res = np.asarray(rouse_sqt_hiq(u))
		return res
		return asarray(rouse_sqt_hiq(u))

		#

pysen/revision.py

+2 −2

Original line number	Diff line number	Diff line
		@@ -3,8 +3,8 @@ PySEN revision module
		"""
		import sys
		__version__ = "2.0"
		__release__ = "a9"
		__date__ = "Dec 19, 2024"
		__release__ = "b1"
		__date__ = "Dec 20, 2024"

		def version(full=False):
		"get pysen version number"

test/test_nscatt.py

+42 −31

Original line number	Diff line number	Diff line
		@@ -37,11 +37,27 @@ class NeutronScatteringTestCase(unittest.TestCase):
		fx = f_debye(inp[:,0])
		self.assertTrue(np.allclose(fx, inp[:,1], atol=1e-03))

		# small argument test
		for x in np.logspace(-13,-3,11):
		self.assertAlmostEqual(1.0, f_debye(x), 5, f"f_Debye for x={x}")

		#
		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_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 = h_rouse_hiq(inp[:,0])
		self.assertTrue(np.allclose(hu, inp[:,1], atol=1e-03))



		def test_rouse_parameters(self):
		"test for Rouse model parameters"
		@@ -65,18 +81,21 @@ class NeutronScatteringTestCase(unittest.TestCase):
		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 = h_rouse_hiq(inp[:,0])
		self.assertTrue(np.allclose(hu, inp[:,1], atol=1e-03))

		def test_rouse_sqt1(self):
		def test_debye_sq(self):
		"test for Debye Sq"
		# Q S(Q)
		d = np.asarray( [[ 0.000, 100.0000],
		[ 0.100, 24.1731],
		[ 0.200, 1.7424],
		[ 0.300, 0.3466],
		[ 0.400, 0.1098],
		[ 0.500, 0.0450]] )
		rmod = RouseModel.from_wl4(Wl4=39500.0, RE=40.0, N=100, T=300)
		for q, expected in d:
		self.assertAlmostEqual(expected, rmod.sq(q), 4, f"SQ(Debye) for q={q}")

		def test_rouse_sqt_datreat(self):
		"test for Rouse Sqt (vs datreat)"
		# results from datreat
		# tau S(Q) S(Q,t)
		@@ -99,9 +118,8 @@ class NeutronScatteringTestCase(unittest.TestCase):
		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"
		def test_rouse_sqt_pep(self):
		"test for Rouse Sqt for PEP"
		# results for PEP@423
		# S(Q) S(Q,t)
		data = [[5.00501524, 5.00501524],
		@@ -116,19 +134,8 @@ class NeutronScatteringTestCase(unittest.TestCase):
		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):
		def test_sqt_rouse_hiq(self):
		"test for Rouse Sqt (high Q)"
		inp = np.array([
		# omegat S(omegat)
		@@ -140,13 +147,17 @@ class NeutronScatteringTestCase(unittest.TestCase):

		u, expected = inp[:,0], inp[:,1]

		# first direct method
		actual, accuracy = rouse_sqt_hiq(inp[:,0], accuracy=True)
		self.assertTrue(np.allclose(actual, expected, atol=1e-03), "sqt_hiq1 (value)")
		self.assertTrue(np.all(abs(accuracy)<1e-5) , "sqt_hiq1 (errors)")

		# second for RouseModel method
		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]
		actual = rmod.sqt_hiq(taus, q)
		self.assertTrue(np.allclose(actual, expected, atol=1e-4), "rouse_hiq2")


		if __name__ == "__main__":
		unittest.main(exit=False, verbosity=2)