diff --git a/examples/test_cheb_lj_2p_xd.py b/examples/test_cheb_lj_2p_xd.py
index a39e07b7e8b3be6b6ac9fea3951533d17ea734e8..6deb8bf833dbe5b10419812d6adf6d7c5c1f2c16 100644
--- a/examples/test_cheb_lj_2p_xd.py
+++ b/examples/test_cheb_lj_2p_xd.py
@@ -6,7 +6,7 @@ from rotagaporp.chebyshev import Chebyshev
# from rotagaporp.symplectic import Symplectic
from rotagaporp.systems import LennardJones, MPPS1D, MPPS3D
from rotagaporp.mpliouville import MPLPQP
-from rotagaporp.liouville import LiouvillePowersQP
+from rotagaporp.liouville import SPLPQP
propagator = Chebyshev
# propagator = Symplectic
@@ -15,7 +15,7 @@ propagator = Chebyshev
params_2p = {'tstart': 0.0, 'tend': 2.0, 'tstep': 0.05, 'nterms': 4,
'liouville': MPLPQP}
params_1p = {'tstart': 0.0, 'tend': 2.0, 'tstep': 0.05, 'nterms': 4,
- 'liouville': LiouvillePowersQP}
+ 'liouville': SPLPQP}
M = [2.0, 2.0]
syspar_3d = {'symbols': ['0x:z', '1x:z'], 'q0': [[0, 0, 0], [2, 0, 0]],
diff --git a/examples/test_cheb_mo_1p_1d_twostep.py b/examples/test_cheb_mo_1p_1d_twostep.py
index 18129790e9375bbfab62e03e9ab9c6a085866e46..afcd77b0b27228f5032fe31ce408b71fbede804d 100644
--- a/examples/test_cheb_mo_1p_1d_twostep.py
+++ b/examples/test_cheb_mo_1p_1d_twostep.py
@@ -4,8 +4,9 @@ import matplotlib.pyplot as plt
from rotagaporp.chebyshev import Chebyshev, TwoStepChebyshev
from rotagaporp.symplectic import Symplectic
from rotagaporp.systems import Morse
+from liouville import SPLPQP
-params = {'tstart': 0., 'tend': 500., 'tstep': 0.05, 'nterms': 4, 'tqdm': True}
+params = {'tstart': 0., 'tend': 500., 'tstep': 0.05, 'nterms': 3, 'tqdm': True}
syspar0 = {'q0': 3., 'p0': 0.}
system = Morse(**syspar0)
@@ -43,8 +44,9 @@ edrift_ch = ch_prop.er
# two-step propagator
params['q_1'] = syspar0['q0']
params['p_1'] = syspar0['p0']
+params['liouville'] = SPLPQP
params['parity'] = 0
-# params['signed'] = False
+params['signed'] = False
ts_prop = TwoStepChebyshev(syst=system, **params)
ts_prop.propagate()
time_ts, q_ts, p_ts = ts_prop.get_trajectory()
diff --git a/examples/test_cheb_mo_2p_xd.py b/examples/test_cheb_mo_2p_xd.py
index e4d32219f0490f4cf6ce682ba602a2736aaadfb6..2edab76b5abd80c2ac5ebf22927f474313efc753 100644
--- a/examples/test_cheb_mo_2p_xd.py
+++ b/examples/test_cheb_mo_2p_xd.py
@@ -6,7 +6,7 @@ from rotagaporp.chebyshev import Chebyshev
# from rotagaporp.symplectic import Symplectic
from rotagaporp.systems import Morse, MPPS1D, MPPS3D
from rotagaporp.mpliouville import MPLPQP
-from rotagaporp.liouville import LiouvillePowersQP
+from rotagaporp.liouville import SPLPQP
propagator = Chebyshev
# propagator = Symplectic
@@ -15,7 +15,7 @@ propagator = Chebyshev
params_2p = {'tstart': 0.0, 'tend': 10., 'tstep': 0.1, 'nterms': 4,
'liouville': MPLPQP}
params_1p = {'tstart': 0.0, 'tend': 10., 'tstep': 0.1, 'nterms': 4,
- 'liouville': LiouvillePowersQP}
+ 'liouville': SPLPQP}
M = [2.0, 2.0]
syspar_3d = {'symbols': ['0x:z', '1x:z'], 'q0': [[1, 0, 0], [4, 0, 0]],
diff --git a/examples/time_1p_nr_cheb.py b/examples/time_1p_nr_cheb.py
index 6d3d24e47d04cfaa0756fd1a62d0321411d54a8a..85d66f8cc40bf9f18ecce9c0f62514001eb5acab 100644
--- a/examples/time_1p_nr_cheb.py
+++ b/examples/time_1p_nr_cheb.py
@@ -2,7 +2,7 @@
import time
import sympy as sp
import numpy as np
-from rotagaporp.liouville import LiouvillePowersNR, LiouvillePowersQP
+from rotagaporp.liouville import SPLPNR, SPLPQP
from rotagaporp.systems import Harmonic, Anharmonic, Morse, LennardJones
nterms = 16
@@ -11,14 +11,14 @@ qval = 5.
pval = -1.
ts = time.time()
-lp_nr = LiouvillePowersNR(nterms, syst, recursion='chebyshev', sign=1)
+lp_nr = SPLPNR(nterms, syst, recursion='chebyshev', sign=1)
te = time.time()
-print('preparation time for LiouvillePowersNR:', te-ts)
+print('preparation time for SPLPNR:', te-ts)
ts = time.time()
-lp_qp = LiouvillePowersQP(nterms, syst, recursion='chebyshev', sign=1)
+lp_qp = SPLPQP(nterms, syst, recursion='chebyshev', sign=1)
te = time.time()
-print("preparation time for LiouvillePowersQP:", te-ts)
+print("preparation time for SPLPQP:", te-ts)
qvec_nr, pvec_nr = lp_nr.get_val_float(qval, pval)
qvec_qp, pvec_qp = lp_qp.get_val_float(qval, pval)
@@ -32,10 +32,10 @@ ts = time.time()
for q, p in zip(qarr, parr):
lp_nr.get_val_float(q, p)
te = time.time()
-print("propagation time for LiouvillePowersNR:", te-ts)
+print("propagation time for SPLPNR:", te-ts)
ts = time.time()
for q, p in zip(qarr, parr):
lp_qp.get_val_float(q, p)
te = time.time()
-print("propagation time for LiouvillePowersQP:", te-ts)
+print("propagation time for SPLPQP:", te-ts)
diff --git a/examples/time_1p_nr_newt.py b/examples/time_1p_nr_newt.py
index 52db0b8db37eea7da96898ef0a38979512f20373..cfca76a7d094b4563f005406c0b93ba666127d28 100644
--- a/examples/time_1p_nr_newt.py
+++ b/examples/time_1p_nr_newt.py
@@ -2,7 +2,7 @@
import time
import sympy as sp
import numpy as np
-from rotagaporp.liouville import LiouvillePowersNR, LiouvillePowersQP
+from rotagaporp.liouville import SPLPNR, SPLPQP
from rotagaporp.systems import Harmonic, Anharmonic, Morse, LennardJones
from rotagaporp.newtonian import leja_points_ga
@@ -14,16 +14,16 @@ pval = -1.
lambdas, scaling = leja_points_ga(nterms)
ts = time.time()
-lp_nr = LiouvillePowersNR(nterms, syst, recursion='newton', nodes=lambdas,
+lp_nr = SPLPNR(nterms, syst, recursion='newton', nodes=lambdas,
scale=1./scaling)
te = time.time()
-print('preparation time for LiouvillePowersNR:', te-ts)
+print('preparation time for SPLPNR:', te-ts)
ts = time.time()
-lp_qp = LiouvillePowersQP(nterms, syst, recursion='newton', nodes=lambdas,
+lp_qp = SPLPQP(nterms, syst, recursion='newton', nodes=lambdas,
scale=1./scaling)
te = time.time()
-print("preparation time for LiouvillePowersQP:", te-ts)
+print("preparation time for SPLPQP:", te-ts)
assert all(nr.equals(qp) for nr, qp in zip(lp_nr.ilnpexpr, lp_qp.pobj.Lf))
@@ -39,10 +39,10 @@ ts = time.time()
for q, p in zip(qarr, parr):
lp_nr.get_val_float(q, p)
te = time.time()
-print("propagation time for LiouvillePowersNR:", te-ts)
+print("propagation time for SPLPNR:", te-ts)
ts = time.time()
for q, p in zip(qarr, parr):
lp_qp.get_val_float(q, p)
te = time.time()
-print("propagation time for LiouvillePowersQP:", te-ts)
+print("propagation time for SPLPQP:", te-ts)
diff --git a/rotagaporp/chebyshev.py b/rotagaporp/chebyshev.py
index eea2d63f41d725f69597c2542a51f43434fdde32..e3779b6c50f94ffc0c7111c67f1b5eeb51f32b45 100644
--- a/rotagaporp/chebyshev.py
+++ b/rotagaporp/chebyshev.py
@@ -101,22 +101,9 @@ class TwoStepChebyshev(TwoStepPolynomialPropagator):
self.nterms = int(np.ceil(alpha)) if nterms is None else nterms
self.efficiency = alpha/self.nterms
bf = besself(self.fpprec, signed, prec=self.prec)
+ assert isinstance(signed, bool)
sign = -1 if signed else 1
assert self.substep is None, 'substep not implemented'
self.coeffs = 2.*cheb_coef(bf, self.nterms, alpha)[self.parity::2]
- self.set_polynomial(system=self.syst, scale=2./self.delta,
- recursion='chebyshev', sign=sign)
-
- def step_nopbc(self, time, q, p):
- """ perform a time step with initial coordinate and momentum """
- # the get_val_float() method has to be adapted to evaluate only
- # either the odd or the even powers
- qvecaux, pvecaux = self.iln.get_val_float(q, p)
- qvec = qvecaux[self.parity::2]
- pvec = pvecaux[self.parity::2]
- #
- qt = np.dot(self.coeffs, qvec) + self.prefact*self.q_1
- pt = np.dot(self.coeffs, pvec) + self.prefact*self.p_1
- self.q_1 = q
- self.p_1 = p
- return (qt, pt)
+ self.set_polynomial(system=self.syst, scale=2./self.delta, sign=sign,
+ recursion='chebyshev', parity=self.parity)
diff --git a/rotagaporp/liouville.py b/rotagaporp/liouville.py
index 5493767d44eed2940513ff5ffc80d0a9720c25f3..78718352019be7a5cb12d54765fd9510562e8b81 100644
--- a/rotagaporp/liouville.py
+++ b/rotagaporp/liouville.py
@@ -1,6 +1,7 @@
""" Class to evaluate the powers of the classical Liouville operator iL """
import sympy as sp
import numpy as np
+from generic import _float
q = sp.Symbol('q')
p = sp.Symbol('p')
@@ -99,6 +100,7 @@ class LiouvillePowers:
class LiouvillePowersQP:
+ """ One-particle Liouville applied to phase space variables q and p """
def __init__(self, nterms, system, **kwargs):
self.qobj = LiouvillePowers(nterms, system.hamiltonian, system.q,
system.p, function=system.q, **kwargs)
@@ -106,11 +108,123 @@ class LiouvillePowersQP:
system.p, function=system.p, **kwargs)
def get_val_float(self, qv, pv):
- """ evaluates (iL)^n q and (iL)^n p for given values of q and p """
+ """ 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 """
-class LiouvillePowersNR(LiouvillePowers):
+ def __init__(self, nterms, hamiltonian, coordinate, momentum, function,
+ scale=1, shift=0, 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
+ * coordinate: coordinate symbol as in the Hamiltonian expression
+ * 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
+ > sign = -1 for Chebyshev polynomials with real argument
+ * parity: parity of the Chebyshev polynomials: 0: even, 1: odd """
+
+ self.nterms = nterms
+ self.H = hamiltonian
+ self.q = coordinate
+ self.p = momentum
+ self.f = function
+ self.fpprec = fpprec
+ self.sign = sign
+ self.scale = scale
+ self.nodes = nodes
+ self.parity = parity
+ self.A = sp.diff(self.H, self.q).doit()*self.scale
+ self.B = sp.diff(self.H, self.p).doit()*self.scale
+ getattr(self, recursion)()
+
+ def poisson(self, phi):
+ """ evaluate the Poisson bracket applied to variable phi, iL.phi """
+ alf = sp.diff(phi, self.q).doit()
+ blf = sp.diff(phi, self.p).doit()
+ return self.B*alf - self.A*blf
+
+ def chebyshev(self):
+ """ evaluate the Chebyshev polynomials with argument iL """
+ assert self.sign in (-1, 1)
+ phi0 = self.f
+ phi1 = self.poisson(self.f)
+ if self.parity is None:
+ self.Lf = [phi0, phi1]
+ for k in range(2, self.nterms):
+ phik = 2*self.poisson(self.Lf[k-1]) + self.sign*self.Lf[k-2]
+ self.Lf.append(phik.expand())
+ else:
+ assert self.parity in (0, 1)
+ assert (self.nterms-1)%2 == self.parity
+ if self.parity == 0:
+ phi2 = 2*self.poisson(phi1) + self.sign*self.f
+ self.Lf = [phi0, phi2]
+ else:
+ 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]
+ self.Lf.append(psik.expand())
+
+ def newton(self):
+ """ evaluate the Newton basis polynomials with argument iL """
+ assert self.nodes is not None
+ self.Lf = [self.f]
+ for n in range(1, self.nterms):
+ expr = self.poisson(self.Lf[n-1]) - self.nodes[n-1]*self.Lf[n-1]
+ self.Lf.append(expr.expand())
+
+ def taylor(self):
+ """ evaluate the monomials (iL)^n.f for the Liouville operator L """
+ self.Lf = [self.f]
+ for n in range(1, self.nterms):
+ self.Lf.append(self.poisson(self.Lf[n-1]).expand())
+
+ def get_val_float(self, qval, pval):
+ """ evaluate the symbolic expressions for given values of q and p """
+ repl = {self.q: qval, self.p: pval}
+ auxfunc = _float(self.fpprec, lambda x, y: x.xreplace(y))
+ return [auxfunc(expr, repl) for expr in self.Lf]
+
+ def __repr__(self):
+ """ string representation of the object in terms of A_k and B_k """
+ A = [sp.diff(self.H, self.q, n) for n in range(self.nterms)]
+ B = [sp.diff(self.H, self.p, n) for n in range(self.nterms)]
+ repl = [(sp.S.Zero, sp.Symbol('0'))]
+ for n in reversed(range(self.nterms)):
+ repl.append((A[n], sp.Symbol('A'+str(n))))
+ repl.append((B[n], sp.Symbol('B'+str(n))))
+ if self.f != self.q and self.f != self.p:
+ aa = [sp.diff(self.f, self.q, n) for n in range(self.nterms)]
+ bb = [sp.diff(self.f, self.p, n) for n in range(self.nterms)]
+ for n in reversed(range(1, self.nterms)):
+ repl.append((aa[n], sp.Symbol('a'+str(n))))
+ repl.append((bb[n], sp.Symbol('b'+str(n))))
+ return 'P_n (f: %s): %s' % (self.f, [e.subs(repl) for e in self.Lf])
+
+
+class SPLPQP:
+ """ Polynomials / monomials of SP Liouville operator applied to q and p """
+ def __init__(self, nterms, system, **kwargs):
+ self.qobj = SPLP(nterms, system.hamiltonian, system.q, system.p,
+ function=system.q, **kwargs)
+ self.pobj = SPLP(nterms, system.hamiltonian, system.q, system.p,
+ function=system.p, **kwargs)
+
+ def get_val_float(self, q, p):
+ """ 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))
+
+
+class SPLPNR(SPLP):
""" use the identity (iL)^n q == (iL)^{n-1} p which is valid if
H(q, p) = p^2/2 + V(q) and non-recursive Chebyshev polynomial
expansion or Newtonian polynomial interpolation """
diff --git a/rotagaporp/propagators.py b/rotagaporp/propagators.py
index b8c2e8b50ccc8a4b7af802333cf6e24982651dda..6d1ad8442cf7edcc1256534c251f62a1eecd634c 100644
--- a/rotagaporp/propagators.py
+++ b/rotagaporp/propagators.py
@@ -90,12 +90,12 @@ class Propagator:
self.er = np.array(func(self.en, self.en0), dtype=float)
-from liouville import LiouvillePowersQP
+from liouville import SPLPQP
class PolynomialPropagator(Propagator):
""" Base class for polynomial classical propagators """
name = 'polynomial'
- def __init__(self, liouville=LiouvillePowersQP, liou_params={},
+ def __init__(self, liouville=SPLPQP, liou_params={},
restart=False, restart_file='propag.pkl', fpprec='fixed',
prec=None, **kwargs):
super().__init__(**kwargs)
@@ -152,9 +152,7 @@ class PolynomialPropagator(Propagator):
class TwoStepPolynomialPropagator(PolynomialPropagator):
""" Base class for two-step polynomial classical propagators """
-
name = 'two-step polynomial'
-
def __init__(self, q_1, p_1, parity=0, **kwargs):
assert parity in (0, 1)
super().__init__(**kwargs)
@@ -162,3 +160,12 @@ class TwoStepPolynomialPropagator(PolynomialPropagator):
self.prefact = 2*parity-1
self.q_1 = q_1
self.p_1 = p_1
+
+ def step_nopbc(self, time, q, p):
+ """ perform a time step with initial coordinate and momentum """
+ qvec, pvec = self.iln.get_val_float(q, p)
+ qt = np.dot(self.coeffs, qvec) + self.prefact*self.q_1
+ pt = np.dot(self.coeffs, pvec) + self.prefact*self.p_1
+ self.q_1 = q
+ self.p_1 = p
+ return (qt, pt)
diff --git a/tests/sp_tests.py b/tests/sp_tests.py
index 9f31b5d0c4d034b9f35eecf854306c147b1d5110..949e03c632ede85000fed18776803f464a78ed3e 100644
--- a/tests/sp_tests.py
+++ b/tests/sp_tests.py
@@ -3,7 +3,8 @@ import unittest
import numpy as np
import sympy as sp
from rotagaporp.liouville import q, p
-from rotagaporp.liouville import LiouvillePowers, LiouvillePowersQP, LiouvillePowersNR
+from rotagaporp.liouville import LiouvillePowers, LiouvillePowersQP
+from rotagaporp.liouville import SPLP, SPLPQP, SPLPNR
from rotagaporp.symplectic import Symplectic
from rotagaporp.chebyshev import Chebyshev, TwoStepChebyshev
from rotagaporp.newtonian import Newtonian
@@ -134,8 +135,8 @@ class LiouvillePowersTest(unittest.TestCase):
assert direct == recurs
-class LiouvillePowersNRTest(unittest.TestCase):
- """ Test the LiouvillePowersNR class """
+class SPLPNRTest(unittest.TestCase):
+ """ Test the SPLPNR class """
def test_nrcheb_symbolic(self):
""" recursive and non-recursive Chebyshev expansions: symbolic """
@@ -145,11 +146,11 @@ class LiouvillePowersNRTest(unittest.TestCase):
p = p
hamiltonian = p**2/2 + sp.Function('V', real=True)(q)
syst = system()
- lp_nr = LiouvillePowersNR(nterms, syst, recursion='chebyshev', sign=1)
- lp_tr = LiouvillePowersQP(nterms, syst, recursion='taylor')
+ lp_nr = SPLPNR(nterms, syst, recursion='chebyshev', sign=1)
+ lp_tr = SPLPQP(nterms, syst, recursion='taylor')
self.assertListEqual(lp_tr.qobj.Lf, lp_nr.ilnq)
self.assertListEqual(lp_tr.pobj.Lf, lp_nr.ilnp)
- lp_cr = LiouvillePowersQP(nterms, syst, recursion='chebyshev', sign=1)
+ lp_cr = SPLPQP(nterms, syst, recursion='chebyshev', sign=1)
self.assertListEqual(lp_cr.qobj.Lf, lp_nr.ilnqexpr.tolist())
self.assertListEqual(lp_cr.pobj.Lf, lp_nr.ilnpexpr.tolist())
@@ -157,8 +158,8 @@ class LiouvillePowersNRTest(unittest.TestCase):
""" recursive and non-recursive Chebyshev expansions: numeric """
nterms = 12
syst = Harmonic()
- lp_nr = LiouvillePowersNR(nterms, syst, recursion='chebyshev', sign=1)
- lp_cr = LiouvillePowersQP(nterms, syst, recursion='chebyshev', sign=1)
+ lp_nr = SPLPNR(nterms, syst, recursion='chebyshev', sign=1)
+ lp_cr = SPLPQP(nterms, syst, recursion='chebyshev', sign=1)
qptest = (1., -1.)
qvec_nr, pvec_nr = lp_nr.get_val_float_direct(*qptest)
qvec_cr, pvec_cr = lp_cr.get_val_float(*qptest)
@@ -177,13 +178,13 @@ class LiouvillePowersNRTest(unittest.TestCase):
hamiltonian = p**2/2 + sp.Function('V', real=True)(q)
syst = system()
lambdas, scaling = leja_points_ga(nterms)
- lp_nr = LiouvillePowersNR(nterms, syst, recursion='newton',
+ lp_nr = SPLPNR(nterms, syst, recursion='newton',
nodes=lambdas, scale=1./scaling)
- lp_tr = LiouvillePowersQP(nterms, syst, recursion='taylor',
+ 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 = LiouvillePowersQP(nterms, syst, recursion='newton',
+ lp_cr = SPLPQP(nterms, syst, recursion='newton',
nodes=lambdas, scale=1./scaling)
qchck = all(n.equals(r) for n, r in zip(lp_nr.ilnqexpr, lp_cr.qobj.Lf))
self.assertTrue(qchck)
@@ -195,9 +196,9 @@ class LiouvillePowersNRTest(unittest.TestCase):
nterms = 12
syst = LennardJones()
lambdas, scaling = leja_points_ga(nterms)
- lp_nr = LiouvillePowersNR(nterms, syst, recursion='newton',
+ lp_nr = SPLPNR(nterms, syst, recursion='newton',
nodes=lambdas, scale=1./scaling)
- lp_cr = LiouvillePowersQP(nterms, syst, recursion='newton',
+ lp_cr = SPLPQP(nterms, syst, recursion='newton',
nodes=lambdas, scale=1./scaling)
qptest = (5., -1.)
qvec_nr, pvec_nr = lp_nr.get_val_float_direct(*qptest)
@@ -209,6 +210,146 @@ class LiouvillePowersNRTest(unittest.TestCase):
self.assertTrue(np.allclose(pvec_nr, pvec_cr))
+class SPLPTest(unittest.TestCase):
+ """ Tests for the SPLP class """
+
+ hamsplit = sp.Function('V', real=True)(q)+sp.Function('T', real=True)(p)
+ hamgener = sp.Function('H', real=True)(q, p)
+ funsplit = sp.Function('fq', real=True)(q)+sp.Function('fp', real=True)(p)
+ fungener = sp.Function('f', real=True)(q, p)
+
+ def test_liouvillepower_expr(self):
+ """ expression tests """
+ ilnp = SPLP(4, hamiltonian=(q**2+p**2)/2, coordinate=q, momentum=p,
+ function=p)
+ self.assertEqual(ilnp.Lf, [p, -q, -p, q])
+ ilnq = SPLP(4, hamiltonian=(q**2+p**2)/2, coordinate=q, momentum=p,
+ function=q)
+ self.assertEqual(ilnq.Lf, [q, p, -q, -p])
+ ilnf = SPLP(4, hamiltonian=(q**2+p**2)/2, coordinate=q, momentum=p,
+ function=p)
+ self.assertEqual(ilnf.Lf, [p, -q, -p, q])
+ ilnf = SPLP(2, hamiltonian=(q**2+p**2)/2, coordinate=q, momentum=p,
+ function=q)
+ self.assertEqual(ilnf.Lf, [q, p])
+
+ def test_liouvillepower_eval(self):
+ """ numeric evaluation tests """
+ ilnp = SPLP(4, hamiltonian=(q**2+p**2)/2, coordinate=q, momentum=p,
+ function=p)
+ self.assertEqual(ilnp.get_val_float(1, 1), [1, -1, -1, 1])
+ hamiltonian = (p**2 + q**4)/2 - q**2
+ ilnq = SPLP(4, hamiltonian=hamiltonian, coordinate=q, momentum=p,
+ function=q)
+ ref_expr = [q, p, -2*q**3 + 2*q, -6*p*q**2 + 2*p]
+ ref_val = [1, 1, 0, -4]
+ self.assertEqual(ilnq.Lf, ref_expr)
+ self.assertEqual(ilnq.get_val_float(1, 1), ref_val)
+
+ def test_liouvillepower_repr_qp(self):
+ """ symbols output tests for iL applied to variables q and p """
+ ilnq = SPLP(6, self.hamsplit, coordinate=q, momentum=p, function=q)
+ ilnp = SPLP(6, self.hamsplit, coordinate=q, momentum=p, function=p)
+ ilnq_ref = (r'P_n (f: q): [q, B1, -A1*B2, A1**2*B3 - A2*B1*B2, '
+ r'-A1**3*B4 + 3*A1*A2*B1*B3 + A1*A2*B2**2 - A3*B1**2*B2, '
+ r'A1**4*B5 - 6*A1**2*A2*B1*B4 - 5*A1**2*A2*B2*B3 + '
+ r'4*A1*A3*B1**2*B3 + 3*A1*A3*B1*B2**2 + 3*A2**2*B1**2*B3 '
+ r'+ A2**2*B1*B2**2 - A4*B1**3*B2]')
+ ilnp_ref = (r'P_n (f: p): [p, -A1, -A2*B1, A1*A2*B2 - A3*B1**2, '
+ r'-A1**2*A2*B3 + 3*A1*A3*B1*B2 + A2**2*B1*B2 - A4*B1**3, '
+ r'A1**3*A2*B4 - 4*A1**2*A3*B1*B3 - 3*A1**2*A3*B2**2 - '
+ r'3*A1*A2**2*B1*B3 - A1*A2**2*B2**2 + 6*A1*A4*B1**2*B2 + '
+ r'5*A2*A3*B1**2*B2 - A5*B1**4]')
+ self.assertEqual(repr(ilnq), ilnq_ref)
+ self.assertEqual(repr(ilnp), ilnp_ref)
+
+ def test_liouvillepower_repr_f(self):
+ """ symbols output tests for iL applied to f = fq(q) + fp(p) """
+ ilnf = SPLP(5, self.hamsplit, coordinate=q, momentum=p,
+ function=self.funsplit)
+ ilnf_ref = (r'P_n (f: fp(p) + fq(q)): [fp(p) + fq(q), -A1*b1 + B1*a1'
+ r', A1**2*b2 - A1*B2*a1 - A2*B1*b1 + B1**2*a2, -A1**3*b3 '
+ r'+ A1**2*B3*a1 + 3*A1*A2*B1*b2 + A1*A2*B2*b1 - 3*A1*B1*B'
+ r'2*a2 - A2*B1*B2*a1 - A3*B1**2*b1 + B1**3*a3, A1**4*b4 -'
+ r' A1**3*B4*a1 - 6*A1**2*A2*B1*b3 - 4*A1**2*A2*B2*b2 - A1'
+ r'**2*A2*B3*b1 + 4*A1**2*B1*B3*a2 + 3*A1**2*B2**2*a2 + 3*'
+ r'A1*A2*B1*B3*a1 + A1*A2*B2**2*a1 + 4*A1*A3*B1**2*b2 + 3*'
+ r'A1*A3*B1*B2*b1 - 6*A1*B1**2*B2*a3 + 3*A2**2*B1**2*b2 + '
+ r'A2**2*B1*B2*b1 - 4*A2*B1**2*B2*a2 - A3*B1**2*B2*a1 - A4'
+ r'*B1**3*b1 + B1**4*a4]')
+ self.assertEqual(repr(ilnf), ilnf_ref)
+
+ def test_liouvillepower_symb(self):
+ """ numeric evaluation tests """
+ ilnq = SPLP(4, hamiltonian=(q**2+p**2)/2, coordinate=q, momentum=p,
+ function=q)
+ self.assertEqual(ilnq.get_val_float(1, 1), [1, 1, -1, -1])
+
+ def test_qp_symmetry(self):
+ """ test iL^n q == iL^(n-1) p when kinetic energy is T(p) = p^2/2 """
+
+ nterms = 5
+ ham = p**2/2 + sp.Function('V', real=True)(q)
+ lpq_t = SPLP(nterms, ham, coordinate=q, momentum=p, function=q,
+ recursion='taylor')
+ lpp_t = SPLP(nterms, ham, coordinate=q, momentum=p, function=p,
+ recursion='taylor')
+ for n in range(1, nterms):
+ self.assertEqual(lpq_t.Lf[n], lpp_t.Lf[n-1])
+
+ # no symmetry if kinetic energy not explicitly defined, T(p)
+ lpq_t = SPLP(nterms, self.hamsplit, coordinate=q, momentum=p,
+ function=q, recursion='taylor')
+ lpp_t = SPLP(nterms, self.hamsplit, coordinate=q, momentum=p,
+ function=p, recursion='taylor')
+ for n in range(1, nterms):
+ self.assertNotEqual(lpq_t.Lf[n], lpp_t.Lf[n-1])
+
+ def test_direct_recursive(self):
+ """ check that the Chebyshev recurrence relations
+ T_n(iL) f = 2 iL T_{n-1}(iL) f -/+ T_{n-2}(iL) f and the expansion
+ T_n(iL) f = sum_{k=0}^n C^(n,-/+)_k (iL)^k f are equal """
+
+ nterms = 5
+ cnk = cnk_chebyshev_T2
+ ilnt = SPLP(nterms, self.hamgener, coordinate=q, momentum=p,
+ function=self.fungener, recursion='taylor')
+ for signed in (True, False):
+ sign = -1 if signed else 1
+ ilnc = SPLP(nterms, self.hamgener, coordinate=q, momentum=p,
+ function=self.fungener, recursion='chebyshev',
+ sign=sign)
+ for n in range(nterms):
+ direct = sum(ilnt.Lf[k]*cnk(n, k, signed) for k in range(n+1))
+ self.assertEqual(direct, ilnc.Lf[n])
+
+ def test_chebyshev_parity(self):
+ """ compare Chebyshev polynomials of same parity with and without
+ parity parameter for different sign, valid for any H and f """
+
+ nterms = 7
+ # 102 seconds for nterms=5
+ # hamiltonian, function = self.hamgener, self.fungener
+ # 11 seconds for nterms=5
+ # hamiltonian, function = self.hamsplit, self.funsplit
+ # 3 seconds for nterms=7
+ hamiltonian, function = self.hamsplit, q
+ for sign in (-1, 1):
+ ilnall = SPLP(nterms, hamiltonian, q, p, function,
+ recursion='chebyshev', sign=sign, parity=None)
+ ilneven = SPLP(nterms, hamiltonian, q, p, function,
+ recursion='chebyshev', sign=sign, parity=0)
+ self.assertEqual(ilnall.Lf[0].factor(), ilneven.Lf[0].factor())
+ self.assertEqual(ilnall.Lf[2].factor(), ilneven.Lf[1].factor())
+ self.assertEqual(ilnall.Lf[4].factor(), ilneven.Lf[2].factor())
+ self.assertEqual(ilnall.Lf[6].factor(), ilneven.Lf[3].factor())
+ ilnodd = SPLP(nterms-1, hamiltonian, q, p, function,
+ recursion='chebyshev', sign=sign, parity=1)
+ self.assertEqual(ilnall.Lf[1].factor(), ilnodd.Lf[0].factor())
+ self.assertEqual(ilnall.Lf[3].factor(), ilnodd.Lf[1].factor())
+ self.assertEqual(ilnall.Lf[5].factor(), ilnodd.Lf[2].factor())
+
+
class ChebyshevPropagatorTest(unittest.TestCase):
""" Tests for the Chebyshev polynomial propagator """
@@ -309,10 +450,10 @@ class ChebyshevPropagatorTest(unittest.TestCase):
self.assertTrue(np.allclose(ptu, pts, rtol=1e-7))
def chebyshev_propagator_nonrecursive_test(self):
- """ test Chebyshev propagator with non-recursive LiouvillePowersNR """
+ """ test Chebyshev propagator with non-recursive SPLPNR """
syst = Harmonic(**self.syspar)
params = self.params.copy()
- params['liouville'] = LiouvillePowersNR
+ params['liouville'] = SPLPNR
prop_nr = Chebyshev(syst=syst, **params)
prop_nr.propagate()
prop_nr.analyse()
@@ -326,7 +467,8 @@ class ChebyshevPropagatorTest(unittest.TestCase):
class TwoStepChebyshevPropagatorTest(unittest.TestCase):
""" Tests for the two-step Chebyshev polynomial propagator """
- params = {'tstart': 0.0, 'tend': 5.0, 'tstep': 0.2, 'nterms': 8}
+ params = {'tstart': 0.0, 'tend': 5.0, 'tstep': 0.2, 'nterms': 7,
+ 'liouville': SPLPQP}
syspar = {'q0': 1.1, 'p0': 0.0}
def chebyshev_propagator_ballistic_and_harmonic_test(self):
@@ -355,10 +497,12 @@ class TwoStepChebyshevPropagatorTest(unittest.TestCase):
params = self.params.copy()
params['q_1'] = qinit[0]
params['p_1'] = pinit[0]
- for parity in (0, 1):
+ for parity, nterms in zip((0, 1), (7, 8)):
for signed in (True, False):
- prop = TwoStepChebyshev(syst=syst, parity=parity,
- signed=signed, **params)
+ params['nterms'] = nterms
+ params['parity'] = parity
+ params['signed'] = signed
+ prop = TwoStepChebyshev(syst=syst, **params)
prop.propagate()
prop.analyse()
self.assertLess(np.max(np.fabs(prop.er)), 1e-8)
@@ -370,7 +514,7 @@ class TwoStepChebyshevPropagatorTest(unittest.TestCase):
def chebyshev_propagator_lj_test(self):
""" test with a particle in Lennard-Jones potential """
syspar = {'p0': -1.0, 'q0': 2}
- params = {'tstart': 0.0, 'tend': 1.8, 'tstep': 0.02, 'nterms': 6}
+ params = {'tstart': 0.0, 'tend': 1.8, 'tstep': 0.02, 'nterms': 5}
syst = LennardJones(**syspar)
ref_prop = Chebyshev(syst=syst, **params)
ref_prop.propagate()
@@ -379,10 +523,13 @@ class TwoStepChebyshevPropagatorTest(unittest.TestCase):
params['q_1'] = q_ref[0]
params['p_1'] = p_ref[0]
params['tstart'] = params['tstart'] + params['tstep']
- for parity in (0, 1):
+ params['liouville'] = SPLPQP
+ for parity, nterms in zip((0, 1), (5, 6)):
for signed in (False, True):
- prop = TwoStepChebyshev(syst=syst, parity=parity,
- signed=signed, **params)
+ params['nterms'] = nterms
+ params['parity'] = parity
+ params['signed'] = signed
+ prop = TwoStepChebyshev(syst=syst, **params)
prop.propagate()
prop.analyse()
self.assertLess(np.max(np.fabs(prop.er)), 1e-3)
@@ -722,8 +869,9 @@ class ChebyshevAssessmentTest(unittest.TestCase):
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,
- 'signed': False, 'tstart': 0., 'tend': 10.0, 'tstep': 0.5}
+ 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()