diff --git a/rotagaporp/chebyshev.py b/rotagaporp/chebyshev.py
index e3779b6c50f94ffc0c7111c67f1b5eeb51f32b45..1179172f38cb15fdf35a6bf24abb80b866c77d64 100644
--- a/rotagaporp/chebyshev.py
+++ b/rotagaporp/chebyshev.py
@@ -56,38 +56,27 @@ def chebyshev_scalar(z, t, delta=None, lambda_min=None, nterms=None,
class Chebyshev(PolynomialPropagator):
""" Chebyshev propagator for classical particle dynamics """
name = 'chebyshev'
- def __init__(self, delta=1., lambda_min=None, nterms=None, signed=False,
- **kwargs):
+ def __init__(self, delta=1., nterms=None, signed=False, **kwargs):
super().__init__(**kwargs)
+ assert delta > 0
self.delta = delta
alpha = self.delta*self.tstep/2.
self.nterms = int(np.ceil(alpha)) if nterms is None else nterms
self.efficiency = alpha/self.nterms
- self.lambda_min = -self.delta/2. if lambda_min is None else lambda_min
- self.scale = 2./self.delta
- self.shift = -(2.*self.lambda_min/self.delta+1.)
bf = besself(self.fpprec, signed, prec=self.prec)
+ assert isinstance(signed, bool)
sign = -1 if signed else 1
- if self.fpprec == 'fixed':
- exp = np.exp
- else:
- exp = lambda x: sp.exp(x).evalf(self.prec)
if self.substep is not None:
- self.cheb_coef = []
+ self.coeffs = []
for stime in self.subtime:
alpha = self.delta*stime/2.
- scaling = exp((self.lambda_min+self.delta/2.0)*stime)
- coef = cheb_coef(bf, self.nterms, alpha)*scaling
- self.cheb_coef.append(coef)
+ self.coeffs.append(cheb_coef(bf, self.nterms, alpha))
else:
- scaling = exp((self.lambda_min+self.delta/2.0)*self.tstep)
- self.cheb_coef = cheb_coef(bf, self.nterms, alpha)*scaling
- self.coeffs = self.cheb_coef
- self.set_polynomial(system=self.syst, scale=self.scale,
- shift=self.shift, recursion='chebyshev',
- sign=sign)
+ self.coeffs = cheb_coef(bf, self.nterms, alpha)
+ self.set_polynomial(system=self.syst, scale=2./self.delta, sign=sign,
+ recursion='chebyshev')
class TwoStepChebyshev(TwoStepPolynomialPropagator):
@@ -96,6 +85,7 @@ class TwoStepChebyshev(TwoStepPolynomialPropagator):
def __init__(self, delta=1., nterms=None, signed=False, **kwargs):
super().__init__(**kwargs)
+ assert delta > 0
self.delta = delta
alpha = self.delta*self.tstep/2.
self.nterms = int(np.ceil(alpha)) if nterms is None else nterms
diff --git a/rotagaporp/liouville.py b/rotagaporp/liouville.py
index 78718352019be7a5cb12d54765fd9510562e8b81..316a23cc645ca8c94e69fd42192210de14c0cf4c 100644
--- a/rotagaporp/liouville.py
+++ b/rotagaporp/liouville.py
@@ -1,4 +1,4 @@
-""" Class to evaluate the powers of the classical Liouville operator iL """
+""" Classes to construct polynomials of the classical Liouville operator """
import sympy as sp
import numpy as np
from generic import _float
@@ -111,11 +111,12 @@ class LiouvillePowersQP:
""" evaluates the expressions for given values of q and p """
return (self.qobj.get_val_float(qv, pv), self.pobj.get_val_float(qv, pv))
+
class SPLP:
""" Polynomials and monomials of single-particle Liouville operator L """
def __init__(self, nterms, hamiltonian, coordinate, momentum, function,
- scale=1, shift=0, recursion='taylor', nodes=None, sign=None,
+ scale=1, recursion='taylor', nodes=None, sign=None,
parity=None, fpprec='fixed'):
""" * nterms: number of polynomials, the highest order (nterms-1)
* hamiltonian: expression of the system Hamiltonian
@@ -123,7 +124,6 @@ class SPLP:
* momentum: momentum symbol as in the Hamiltonian expression
* function: a system variable to which iL is applied
* scale: a scaling factor for L
- * shift: a shift for L that is always 0, deprecated argument
* nodes: a list of Newton interpolation points
* sign: the sign in the Chebyshev polynomial recursion formula
> sign = +1 for Chebyshev polynomials with imaginary argument
@@ -170,8 +170,8 @@ class SPLP:
phi3 = 4*self.poisson(self.poisson(phi1)) + self.sign*3*phi1
self.Lf = [phi1, phi3]
for k in range(2, (self.nterms-1)//2+1):
- Lsq_psi_1 = self.poisson(self.poisson(self.Lf[k-1]))
- psik = 4*Lsq_psi_1 + self.sign*2*self.Lf[k-1] - self.Lf[k-2]
+ sq_psi = self.poisson(self.poisson(self.Lf[k-1]))
+ psik = 4*sq_psi + self.sign*2*self.Lf[k-1] - self.Lf[k-2]
self.Lf.append(psik.expand())
def newton(self):
@@ -219,9 +219,11 @@ class SPLPQP:
self.pobj = SPLP(nterms, system.hamiltonian, system.q, system.p,
function=system.p, **kwargs)
- def get_val_float(self, q, p):
+ def get_val_float(self, qval, pval):
""" evaluate the symbolic expressions for given values of q and p """
- return (self.qobj.get_val_float(q, p), self.pobj.get_val_float(q, p))
+ qvec = self.qobj.get_val_float(qval, pval)
+ pvec = self.pobj.get_val_float(qval, pval)
+ return (qvec, pvec)
class SPLPNR(SPLP):
@@ -232,12 +234,9 @@ class SPLPNR(SPLP):
def __init__(self, nterms, system, recursion=None, **kwargs):
""" the kinetic energy must be explicit and H(q, p) = p^2/2 + V(q) """
from combinatorics import cnk_chebyshev_T2, cnk_newton_1
- from generic import _float
assert recursion in ['chebyshev', 'newton']
self.scale = kwargs.get('scale', 1)
- self.shift = kwargs.get('shift', 0)
- assert np.isclose(self.shift, 0.0), 'currently only without shift'
fpprec = kwargs.get('fpprec', 'fixed')
prec = kwargs.get('prec', None)
# current implementation only for float
@@ -271,10 +270,11 @@ class SPLPNR(SPLP):
self.get_val_float = self.get_val_float_nrnewt
def expand_newt_expr(self, iln):
+ """ evaluate the Newton polynomial using the C^n_k coefficients """
newt = []
- for n in range(len(iln)):
- iln_rev = reversed(iln[:n+1])
- newt.append(sum(self.cnk[n, k]*t for k, t in enumerate(iln_rev)))
+ for num in range(len(iln)):
+ ilnrev = reversed(iln[:num+1])
+ newt.append(sum(self.cnk[num, k]*t for k, t in enumerate(ilnrev)))
return np.array(newt)
def get_val_float_direct(self, qval, pval):
diff --git a/rotagaporp/newtonian.py b/rotagaporp/newtonian.py
index e18adba207475295d07618c7cbf7ba0774053688..fec3389de0c0880d23df5db9403f7393a5370344 100644
--- a/rotagaporp/newtonian.py
+++ b/rotagaporp/newtonian.py
@@ -41,7 +41,7 @@ def trial_points(number, dtype=float, fpprec='fixed', prec=None):
cosf = {'fixed': np.cos, 'multi': lambda x: sp.N(sp.cos(x), prec)}
sinf = {'fixed': np.sin, 'multi': lambda x: sp.N(sp.sin(x), prec)}
fspace = {
- 'fixed': lambda l, h, n: np.linspace(l, h, n),
+ 'fixed': np.linspace,
'multi': lambda l, h, n: np.arange(l, h+(h-l)/(n-1), (h-l)/(n-1))
}
center = zero[fpprec]
@@ -142,35 +142,33 @@ def newton_scalar(lambdas, func, vals, scale=1.0, shift=0.0):
class Newtonian(PolynomialPropagator):
""" Newtonian propagator for classical particle dynamics """
name = 'newtonian'
- def __init__(self, delta=1., lambda_min=None, nterms=None, **kwargs):
+ def __init__(self, delta=1., nterms=None, **kwargs):
super().__init__(**kwargs)
assert delta > 0
- assert nterms is not None and nterms > 0
self.delta = delta
+ assert nterms is not None and nterms > 0
self.nterms = nterms
- self.lambda_min = -self.delta/2. if lambda_min is None else lambda_min
lambdas, scaling = self.ip_points(algo='ga')
- self.scale = self.delta*scaling
- self.shift = (self.lambda_min+self.delta/2.)/self.scale
- lamb_sc = (lambdas*self.scale+self.shift)
+ scale = self.delta*scaling
+ lamb_sc = lambdas*scale
if self.substep is None:
- self.ddiffs = self.dd_coeffs(lambdas, lamb_sc, self.tstep)
+ ddiffs = self.dd_coeffs(lambdas, lamb_sc, self.tstep)
else:
- self.ddiffs = np.array([self.dd_coeffs(lambdas, lamb_sc, st)
- for st in self.subtime])
- self.coeffs = self.ddiffs
- self.set_polynomial(system=self.syst, scale=1./self.scale,
- shift=-self.shift, recursion='newton',
- nodes=lambdas)
+ ddiffs = np.array([self.dd_coeffs(lambdas, lamb_sc, st)
+ for st in self.subtime])
+ self.coeffs = ddiffs
+ self.set_polynomial(system=self.syst, scale=1./scale, nodes=lambdas,
+ recursion='newton')
def dd_coeffs(self, lamb, lamb_sc, t, algo=1):
""" select the algorithm for divided differences """
algo_map = {1: divdiff_1, 2: divdiff_2}
if self.fpprec == 'fixed':
- return algo_map[algo](lamb, np.exp(lamb_sc*t))
+ retval = algo_map[algo](lamb, np.exp(lamb_sc*t))
elif self.fpprec == 'multi':
- return algo_map[algo](lamb, [sp.exp(i) for i in lamb_sc*t])
+ retval = algo_map[algo](lamb, [sp.exp(i) for i in lamb_sc*t])
+ return retval
def ip_points(self, dtype=float, algo='ga'):
""" select the algorithm for interpolation points """
diff --git a/rotagaporp/propagators.py b/rotagaporp/propagators.py
index 6d1ad8442cf7edcc1256534c251f62a1eecd634c..3b0f30d3ad0b15e4c77767844ea2f7c72c6f4a06 100644
--- a/rotagaporp/propagators.py
+++ b/rotagaporp/propagators.py
@@ -1,6 +1,7 @@
""" Propagators for classical particle dynamics """
import numpy as np
+from liouville import SPLPQP
class Propagator:
""" Base class for all propagators for classical particle dynamics """
@@ -90,14 +91,17 @@ class Propagator:
self.er = np.array(func(self.en, self.en0), dtype=float)
-from liouville import SPLPQP
-
class PolynomialPropagator(Propagator):
""" Base class for polynomial classical propagators """
+
name = 'polynomial'
- def __init__(self, liouville=SPLPQP, liou_params={},
- restart=False, restart_file='propag.pkl', fpprec='fixed',
- prec=None, **kwargs):
+ coeffs = None
+ nterms = None
+ iln = None
+
+ def __init__(self, liouville=SPLPQP, liou_params={}, restart=False,
+ restart_file='propag.pkl', fpprec='fixed', prec=None,
+ **kwargs):
super().__init__(**kwargs)
self.step = self.step_pbc if self.pbc else self.step_nopbc
self.liouville = liouville
@@ -134,8 +138,8 @@ class PolynomialPropagator(Propagator):
import cloudpickle as pkl
with open(filename, 'rb') as fh:
that = pkl.load(fh)
- params = ['__class__', 'name', 'nterms', 'delta', 'tstep', 'scale',
- 'shift', 'liouville', 'liou_params', 'fpprec', 'prec']
+ params = ['__class__', 'name', 'nterms', 'delta', 'tstep', 'liouville',
+ 'liou_params', 'fpprec', 'prec']
for key in params:
attr = getattr(that, key)
errstr = key + ': ' + str(getattr(self, key)) + ' != ' + str(attr)
@@ -152,7 +156,11 @@ class PolynomialPropagator(Propagator):
class TwoStepPolynomialPropagator(PolynomialPropagator):
""" Base class for two-step polynomial classical propagators """
+
name = 'two-step polynomial'
+ coeffs = None
+ iln = None
+
def __init__(self, q_1, p_1, parity=0, **kwargs):
assert parity in (0, 1)
super().__init__(**kwargs)
diff --git a/tests/sp_tests.py b/tests/sp_tests.py
index 949e03c632ede85000fed18776803f464a78ed3e..9793ccab99684e32d07f1536938175821789d4e4 100644
--- a/tests/sp_tests.py
+++ b/tests/sp_tests.py
@@ -3,7 +3,7 @@ import unittest
import numpy as np
import sympy as sp
from rotagaporp.liouville import q, p
-from rotagaporp.liouville import LiouvillePowers, LiouvillePowersQP
+from rotagaporp.liouville import LiouvillePowers
from rotagaporp.liouville import SPLP, SPLPQP, SPLPNR
from rotagaporp.symplectic import Symplectic
from rotagaporp.chebyshev import Chebyshev, TwoStepChebyshev
@@ -178,14 +178,13 @@ class SPLPNRTest(unittest.TestCase):
hamiltonian = p**2/2 + sp.Function('V', real=True)(q)
syst = system()
lambdas, scaling = leja_points_ga(nterms)
- lp_nr = SPLPNR(nterms, syst, recursion='newton',
- nodes=lambdas, scale=1./scaling)
- lp_tr = SPLPQP(nterms, syst, recursion='taylor',
- scale=1./scaling)
+ lp_nr = SPLPNR(nterms, syst, recursion='newton', scale=1./scaling,
+ nodes=lambdas)
+ lp_tr = SPLPQP(nterms, syst, recursion='taylor', scale=1./scaling)
self.assertListEqual(lp_tr.qobj.Lf, lp_nr.ilnq)
self.assertListEqual(lp_tr.pobj.Lf, lp_nr.ilnp)
- lp_cr = SPLPQP(nterms, syst, recursion='newton',
- nodes=lambdas, scale=1./scaling)
+ lp_cr = SPLPQP(nterms, syst, recursion='newton', scale=1./scaling,
+ nodes=lambdas)
qchck = all(n.equals(r) for n, r in zip(lp_nr.ilnqexpr, lp_cr.qobj.Lf))
self.assertTrue(qchck)
pchck = all(n.equals(r) for n, r in zip(lp_nr.ilnpexpr, lp_cr.pobj.Lf))
@@ -196,10 +195,10 @@ class SPLPNRTest(unittest.TestCase):
nterms = 12
syst = LennardJones()
lambdas, scaling = leja_points_ga(nterms)
- lp_nr = SPLPNR(nterms, syst, recursion='newton',
- nodes=lambdas, scale=1./scaling)
- lp_cr = SPLPQP(nterms, syst, recursion='newton',
- nodes=lambdas, scale=1./scaling)
+ lp_nr = SPLPNR(nterms, syst, recursion='newton', scale=1./scaling,
+ nodes=lambdas)
+ lp_cr = SPLPQP(nterms, syst, recursion='newton', scale=1./scaling,
+ nodes=lambdas)
qptest = (5., -1.)
qvec_nr, pvec_nr = lp_nr.get_val_float_direct(*qptest)
qvec_cr, pvec_cr = lp_cr.get_val_float(*qptest)
@@ -866,20 +865,6 @@ class ChebyshevAssessmentTest(unittest.TestCase):
self.assertTrue(np.allclose(pt, ptr, atol=1e-3))
self.assertLess(np.max(np.fabs(prop.er)), 1e-3)
- def chebyshev_shift_test(self):
- """ chebyshev propagator with scale and shift """
- syst = Harmonic(q0=1., p0=0.)
- params = {'delta': 3., 'lambda_min': -2., 'nterms': 16, 'tstart': 0.,
- 'tend': 10.0, 'tstep': 0.5, 'liouville': LiouvillePowersQP,
- 'signed': False}
- prop = Chebyshev(syst=syst, **params)
- prop.propagate()
- prop.analyse()
- t, q, p = prop.get_trajectory()
- qr, pr = syst.solution(times=t)
- self.assertTrue(np.allclose(q, qr))
- self.assertLess(np.max(np.fabs(prop.er)), 1e-6)
-
def chebyshev_symplectic_2p_lj_test(self):
""" two particles in Lennard-Jones potential """