diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..c0a5c2ec9b3373353ba072982f10d61e4ac95713
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+TODO
+pylintrc
+ideas
+measurements*
+.venv*
diff --git a/classical/plot_iln_np_3d.py b/classical/plot_iln_np_3d.py
deleted file mode 100644
index 171c95580e89bd405387f2a197ce45eb46da95d2..0000000000000000000000000000000000000000
--- a/classical/plot_iln_np_3d.py
+++ /dev/null
@@ -1,24 +0,0 @@
-from pandasplot import PandasPlot
-
-var_labels = {
-    'nterms': r'$N$',
-    'totals.operations': 'number of operations',
-    'totals.storage': 'number of atomic objects',
-    'preparation time': 'preparation time',
-    'class': 'MPLP class',
-    'nparticles': 'N particles',
-    'variant': 'variant'
-}
-
-pplt = PandasPlot.from_file('exprstats.json', showy='preparation time',
-                            showx='nparticles', labels=var_labels)
-
-legendvars = ['class']
-# legendvars = ['variant']
-pplt.add_datasets(
-    select={'variant': 'default', 'nterms': 5},
-    vary={v: sorted(list(set(pplt.df[v].tolist()))) for v in legendvars}
-)
-
-pplt.set_axis_type(ytype='log')
-pplt.show(legend=True)
diff --git a/classical/time_uniq.py b/classical/time_uniq.py
deleted file mode 100644
index 8df428d81dd6589256d763648f0bc64dc8cbeb9b..0000000000000000000000000000000000000000
--- a/classical/time_uniq.py
+++ /dev/null
@@ -1,16 +0,0 @@
-import numpy as np
-from timing import timeit
-from generic import uniq1, uniq2
-
-@timeit
-def time_uniq1(seq):
-    return uniq1(seq)
-
-@timeit
-def time_uniq2(seq):
-    return uniq2(seq)
-
-for irange in [10, 100]:
-    print('integers range', irange)
-    seq = np.random.randint(0, irange, (1000000, ))
-    assert time_uniq1(seq) == time_uniq2(seq)
diff --git a/docs/.gitignore b/docs/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/examples/requirements.txt b/examples/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7aff1cb889bce8726e823534e5e41572258acd6e
--- /dev/null
+++ b/examples/requirements.txt
@@ -0,0 +1,5 @@
+pandas==0.24.2
+pycairo==1.18.0
+PyQt5==5.12.1
+PyYAML==5.1
+tensorflow==1.13.1
diff --git a/classical/test_cheb_00_2p_1d_pbc.py b/examples/test_cheb_00_2p_1d_pbc.py
similarity index 89%
rename from classical/test_cheb_00_2p_1d_pbc.py
rename to examples/test_cheb_00_2p_1d_pbc.py
index b9881e99f8485d82ae7bccbfd9456a55550684d3..e86ed4e925137b2acb14b7402a343169a42d712e 100644
--- a/classical/test_cheb_00_2p_1d_pbc.py
+++ b/examples/test_cheb_00_2p_1d_pbc.py
@@ -1,9 +1,9 @@
 """ One dimensional ideal gas (with no interactions) with two particles """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import MPPS1D
-from mpliouville import MPLPDV
 import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import MPPS1D
+from rotagaporp.mpliouville import MPLPDV
 
 syspar = {'symbols': ['0', '1'], 'q0': [[1], [3]], 'p0': [[1], [-1]],
           'masses': [1.0, 1.0], 'pbc': True, 'cell': [5.0],
diff --git a/classical/test_cheb_00_2p_3d_pbc.py b/examples/test_cheb_00_2p_3d_pbc.py
similarity index 87%
rename from classical/test_cheb_00_2p_3d_pbc.py
rename to examples/test_cheb_00_2p_3d_pbc.py
index 919562f5131150ebb8a8b72f636cddae6efd9050..1a4ebd2efc61d4db2eb0a36b2bd3ff1dc42c88ec 100644
--- a/classical/test_cheb_00_2p_3d_pbc.py
+++ b/examples/test_cheb_00_2p_3d_pbc.py
@@ -1,12 +1,12 @@
 """ Three dimensional ideal gas (no interactions) with two particles """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import ManyAtomsSystem
-from mpliouville import MPLPDV
-import matplotlib.pyplot as plt
 import numpy as np
 from ase import Atoms
-from units import LJUnitsElement
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from rotagaporp.mpliouville import MPLPDV
+from utilities.units import LJUnitsElement
 
 argon_u = LJUnitsElement('Ar')
 symbols = 'ArAr'
diff --git a/examples/test_cheb_argon.py b/examples/test_cheb_argon.py
new file mode 100644
index 0000000000000000000000000000000000000000..8ca1e6b7ce601bc09643459d6f72d235471a1c43
--- /dev/null
+++ b/examples/test_cheb_argon.py
@@ -0,0 +1,41 @@
+""" Molecular dynamics of argon at constant energy """
+import time
+from ase import units
+from ase.build import bulk
+from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
+from rotagaporp.systems import ManyAtomsSystem
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.mpliouville import MPLPDV
+from utilities.trajectory import write_ase_trajectory
+
+T0_lj = 3.0 # initial temperature in Ar LJ units
+rho_lj = 0.1 # 0.88 # density in Ar LJ units
+sigma_a = 3.405*units.Angstrom # in Angstrom for Ar
+epsilon_a = 1.65e-21*units.J # in Joule for Ar
+T = T0_lj * epsilon_a / units.kB # initial temperature in Kelvin
+rho = rho_lj/sigma_a**3 # density in atoms per cubic Angstrom
+n_cell = 2 # number of unit cells in one direction
+n_atoms = 4*n_cell**3 # total number of atoms, 4 atoms per unit cell
+l_box = (n_atoms/rho)**(1./3) # edge length of simulation box
+l_cell = l_box / n_cell # latice parameter
+
+atoms = bulk('Ar', 'fcc', a=l_cell, cubic=True).repeat((n_cell, n_cell, n_cell))
+rho_real = sum(atoms.get_masses())/atoms.get_volume()*units.m**3/units.kg
+print('density, kg/m^3', rho_real)
+
+prep_start = time.time()
+MaxwellBoltzmannDistribution(atoms, T*units.kB)
+syst = ManyAtomsSystem(atoms, 'FFLennardJones', numeric={'tool': 'theano'})
+params = {'tstart': 0.0, 'tend': 10.0, 'tstep': 0.001, 'tqdm': True,
+          'nterms': 4, 'liouville': MPLPVR, 'pbc': True}
+prop = Chebyshev(syst=syst, **params)
+prep_end = time.time()
+
+print('preparation time:', prep_end-prep_start)
+print('propagation begins')
+prop_start = time.time()
+prop.propagate()
+prop_end = time.time()
+print('propagation time:', prop_end-prop_start)
+prop.analyse(step=10)
+write_ase_trajectory(prop, atoms, step=10)
diff --git a/classical/test_cheb_ho_2p_1d.py b/examples/test_cheb_ho_2p_1d.py
similarity index 86%
rename from classical/test_cheb_ho_2p_1d.py
rename to examples/test_cheb_ho_2p_1d.py
index 95e8625f540d065b3260d53345b42e4a2be72a39..93d51f55cba81e4191d16e737bf9ae8046899287 100644
--- a/classical/test_cheb_ho_2p_1d.py
+++ b/examples/test_cheb_ho_2p_1d.py
@@ -1,25 +1,24 @@
 """ Two particles interacting in harmonic potential """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import Harmonic, MPPS1D
-from mpliouville import MPLPDV
-import matplotlib.pyplot as plt
 import sympy as sp
 import numpy as np
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import Harmonic, MPPS1D
+from rotagaporp.mpliouville import MPLPDV
 
 M = np.full(2, 2.0)
 r0 = 3.5
 kappa = sp.Rational(1, 2)
 syspar = {'symbols': ['0', '1'], 'q0': [[0], [1]], 'p0': [[1], [0]],
-          'masses': M, 'numeric': {'tool': 'theano'}}
+          'masses': M}
 ffpar = [{'class': 'FFHarmonic', 'pair': ((0,), (1,)),
           'params': {'distance': r0, 'kappa': kappa}}]
 syst = MPPS1D(ffpar, **syspar)
 
-params = {'tstart': 0.0, 'tend': 4.0, 'tstep': 0.001, 'nterms': 4, 'liouville': MPLPDV}
+params = {'tstart': 0.0, 'tend': 4.0, 'tstep': 0.05, 'nterms': 4,
+          'liouville': MPLPDV}
 prop = Chebyshev(syst=syst, **params)
-prop.dump()
-exit()
 prop.propagate()
 prop.analyse()
 (ti, qt, pt) = prop.get_trajectory()
@@ -43,7 +42,7 @@ sys_1p = Harmonic(**sys_1p_par)
 qrta = np.real(sys_1p.solution(sys_1p_par['q0'], sys_1p_par['p0'], tir)[0])
 prta = np.real(sys_1p.solution(sys_1p_par['q0'], sys_1p_par['p0'], tir)[1])
 
-# 2p -> 1p reduced chebyshev solution 
+# 2p -> 1p reduced chebyshev solution
 qrtr = np.abs(qt[:, 0]-qt[:, 1]) - r0
 prtr = (pt[:, 1]*M[0]-pt[:, 0]*M[1])/np.sum(M)
 
diff --git a/examples/test_cheb_ho_2p_1d_restart.py b/examples/test_cheb_ho_2p_1d_restart.py
new file mode 100644
index 0000000000000000000000000000000000000000..3d1be862cc9da2260392570d40eb59c7cb476148
--- /dev/null
+++ b/examples/test_cheb_ho_2p_1d_restart.py
@@ -0,0 +1,77 @@
+""" Two particles interacting in harmonic potential """
+import sympy as sp
+import numpy as np
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import Harmonic, MPPS1D
+from rotagaporp.mpliouville import MPLPDV
+
+M = np.full(2, 2.0)
+r0 = 3.5
+kappa = sp.Rational(1, 2)
+syspar = {'symbols': ['0', '1'], 'q0': [[0], [1]], 'p0': [[1], [0]],
+          'masses': M, 'numeric': {'tool': 'subs'}}
+ffpar = [{'class': 'FFHarmonic', 'pair': ((0,), (1,)),
+          'params': {'distance': r0, 'kappa': kappa}}]
+syst = MPPS1D(ffpar, **syspar)
+
+params = {'tstart': 0.0, 'tend': 4.0, 'tstep': 0.05, 'nterms': 4, 'liouville': MPLPDV, 'tqdm': True}
+
+prop1 = Chebyshev(syst=syst, **params)
+prop1.dump()
+
+prop2 = Chebyshev(syst=syst, restart=True, **params)
+prop2.propagate()
+prop2.analyse()
+(ti, qt, pt) = prop2.get_trajectory()
+
+params = {'tstart': 0.0, 'tend': 4.0, 'tstep': 0.05}
+prop_ref = Symplectic(syst=syst, **params)
+prop_ref.propagate()
+prop_ref.analyse()
+(tir, qtr, ptr) = prop_ref.get_trajectory()
+
+# analytic solution of equivalent 1p problem
+sys_1p_par = {
+    'q0': np.abs(syspar['q0'][0][0]-syspar['q0'][1][0])-r0,
+    # from the motion of center of mass
+    'p0': (M[0]*syspar['p0'][1][0]+M[1]*syspar['p0'][0][0])/np.sum(M),
+    # from total kinetic energy = sum of particles's kinetic energies
+    # 'p0': np.sqrt((M[0]*syspar['p0'][1]**2+M[1]*syspar['p0'][0]**2)/np.sum(M)),
+    'mass': np.prod(M)/np.sum(M)
+}
+sys_1p = Harmonic(**sys_1p_par)
+qrta = np.real(sys_1p.solution(sys_1p_par['q0'], sys_1p_par['p0'], tir)[0])
+prta = np.real(sys_1p.solution(sys_1p_par['q0'], sys_1p_par['p0'], tir)[1])
+
+# 2p -> 1p reduced chebyshev solution 
+qrtr = np.abs(qt[:, 0]-qt[:, 1]) - r0
+prtr = (pt[:, 1]*M[0]-pt[:, 0]*M[1])/np.sum(M)
+
+fig = plt.figure()
+plot = fig.add_subplot(311)
+plot.plot(tir, qrta, label='r, analytic 1p')
+plot.plot(tir, qrtr, label='q_r chebyshev 2p->1p reduced')
+plot.plot(tir, prta, label='p_r analytic 1p')
+plot.plot(tir, prtr, label='p_r chebyshev 2p->1p reduced')
+plot.set_ylabel(r'$r,\quad p_r$')
+plt.legend()
+
+plot = fig.add_subplot(312)
+plt.plot(tir, qtr[:, 0], label='q0 symplectic')
+plt.plot(ti, qt[:, 0], label='q0 chebyshev')
+plt.plot(tir, qtr[:, 1], label='q1 symplectic')
+plt.plot(ti, qt[:, 1], label='q1 chebyshev')
+plot.set_ylabel(r'$q_0,\quad q_1$')
+plt.legend()
+
+plot = fig.add_subplot(313)
+plt.plot(tir, ptr[:, 0], label='p0 symplectic')
+plt.plot(ti, pt[:, 0], label='p0 chebyshev')
+plt.plot(tir, ptr[:, 1], label='p1 symplectic')
+plt.plot(ti, pt[:, 1], label='p1 chebyshev')
+plot.set_xlabel('time')
+plot.set_ylabel(r'$p_0,\quad p_1$')
+plt.legend()
+plt.show()
diff --git a/classical/test_cheb_ho_2p_3d.py b/examples/test_cheb_ho_2p_3d.py
similarity index 88%
rename from classical/test_cheb_ho_2p_3d.py
rename to examples/test_cheb_ho_2p_3d.py
index f499f2175fcfd8ab057c64746811a48a69ec31a6..9b82806c0d0b34409c16e056c1b46dcdc1dced96 100644
--- a/classical/test_cheb_ho_2p_3d.py
+++ b/examples/test_cheb_ho_2p_3d.py
@@ -1,11 +1,11 @@
 """ Two particles interacting in harmonic potential """
-from chebyshev import Chebyshev
-from systems import Harmonic, MPPS1D, MPPS3D
-from mpliouville import MPLPQP
-import matplotlib.pyplot as plt
-from ase.geometry import get_distances
 import sympy as sp
 import numpy as np
+import matplotlib.pyplot as plt
+from ase.geometry import get_distances
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.systems import Harmonic, MPPS1D, MPPS3D
+from rotagaporp.mpliouville import MPLPQP
 
 params = {'tstart': 0.0, 'tend': 1.0, 'tstep': 0.1, 'nterms': 4,
           'liouville': MPLPQP, 'tqdm': True}
@@ -23,8 +23,9 @@ prop_1d.propagate()
 prop_1d.analyse()
 (ti1d, qt1d, pt1d) = prop_1d.get_trajectory()
 
-syspar_3d = {'symbols': ['0x:z', '1x:z'], 'q0': [[1, 2, 3], [1.7071, 2.7071, 3]],
-          'p0': [[0.7071, 0.7071, 0], [0, 0, 0]], 'masses': M}
+syspar_3d = {'symbols': ['0x:z', '1x:z'],
+             'q0': [[1, 2, 3], [1.7071, 2.7071, 3]],
+             'p0': [[0.7071, 0.7071, 0], [0, 0, 0]], 'masses': M}
 
 ffpar_3d = [{'class': 'FFHarmonic', 'pair': ((0, 1, 2), (3, 4, 5)),
              'params': {'distance': r0, 'kappa': kappa}}]
diff --git a/classical/test_cheb_lj_2p_1d.py b/examples/test_cheb_lj_2p_1d.py
similarity index 88%
rename from classical/test_cheb_lj_2p_1d.py
rename to examples/test_cheb_lj_2p_1d.py
index 9901574f2f994c20d2cc19651208d38283114527..05518e270d6f2555a07cdecefdcc9d1792e3cb03 100644
--- a/classical/test_cheb_lj_2p_1d.py
+++ b/examples/test_cheb_lj_2p_1d.py
@@ -1,9 +1,9 @@
 """ One dimensional Lennard-Jones gas with two particles """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import MPPS1D
-from mpliouville import MPLPQP
 import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import MPPS1D
+from rotagaporp.mpliouville import MPLPQP
 
 syspar = {'symbols': ['0', '1'], 'q0': [[0], [2]], 'p0': [[0], [-1]],
           'masses': [1.0, 1.0]}
diff --git a/classical/test_cheb_lj_2p_1d_pbc.py b/examples/test_cheb_lj_2p_1d_pbc.py
similarity index 91%
rename from classical/test_cheb_lj_2p_1d_pbc.py
rename to examples/test_cheb_lj_2p_1d_pbc.py
index da7d65b8cbd57665d2e9bd2b74549da56f7982c8..56f1f1db5cfd766716d2d12e4d1eb2d3cb1f56f7 100644
--- a/classical/test_cheb_lj_2p_1d_pbc.py
+++ b/examples/test_cheb_lj_2p_1d_pbc.py
@@ -1,9 +1,9 @@
 """ One dimensional Lennard-Jones gas with two particles """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import MPPS1D
-from mpliouville import MPLPDV
 import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import MPPS1D
+from rotagaporp.mpliouville import MPLPDV
 
 syspar = {'symbols': ['0', '1'], 'q0': [[1], [3]], 'p0': [[0], [-1]],
           'masses': [1.0, 1.0], 'pbc': True, 'cell': [15.],
diff --git a/classical/test_cheb_lj_2p_3d.py b/examples/test_cheb_lj_2p_3d.py
similarity index 85%
rename from classical/test_cheb_lj_2p_3d.py
rename to examples/test_cheb_lj_2p_3d.py
index 06b814ab1313dd1bb62f7e97dd9770a915a440a7..06c85c43be9f8aad501bf8fb2dc25e728bc434dd 100644
--- a/classical/test_cheb_lj_2p_3d.py
+++ b/examples/test_cheb_lj_2p_3d.py
@@ -1,13 +1,13 @@
 """ Three dimensional Lennard-Jones gas with two particles """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import ManyAtomsSystem
-from mpliouville import MPLPNR
-import matplotlib.pyplot as plt
 import numpy as np
 from ase import Atoms
-from units import LJUnitsElement
-from trajectory import write_ase_trajectory
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from rotagaporp.mpliouville import MPLPNR
+from utilities.units import LJUnitsElement
+from utilities.trajectory import write_ase_trajectory
 
 argon_u = LJUnitsElement('Ar')
 symbols = 'ArAr'
@@ -27,7 +27,7 @@ numeric = {'tool': 'lambdify', 'modules': 'tensorflow'}
 # numeric = {'tool': 'subs'}
 syst = ManyAtomsSystem(atoms, 'FFLennardJones', numeric=numeric)
 
-params = {'tstart': 0.0, 'tend': 10.0, 'tstep': 0.001, 'nterms': 4,
+params = {'tstart': 0.0, 'tend': 10.0, 'tstep': 0.005, 'nterms': 4,
           'liouville': MPLPNR, 'tqdm': True}
 prop = Chebyshev(syst=syst, **params)
 print('propagation begins')
@@ -36,7 +36,7 @@ prop.propagate()
 (ti, qt, pt) = prop.get_trajectory_3d()
 write_ase_trajectory(prop, atoms)
 
-params = {'tstart': 0.0, 'tend': 14.0, 'tstep': 0.005, 'tqdm': True}
+params = {'tstart': 0.0, 'tend': 10.0, 'tstep': 0.005, 'tqdm': True}
 prop_ref = Symplectic(syst=syst, **params)
 prop_ref.propagate()
 # prop_ref.analyse()
diff --git a/classical/test_cheb_lj_2p_3d_pbc.py b/examples/test_cheb_lj_2p_3d_pbc.py
similarity index 78%
rename from classical/test_cheb_lj_2p_3d_pbc.py
rename to examples/test_cheb_lj_2p_3d_pbc.py
index 97d7bbc9e652a5cabceab84a5bbdb5bd04317fd4..0cf2d512bb9fd4d28c95e55588d320f18260b6f2 100644
--- a/classical/test_cheb_lj_2p_3d_pbc.py
+++ b/examples/test_cheb_lj_2p_3d_pbc.py
@@ -1,12 +1,12 @@
 """ Three dimensional Lennard-Jones gas with two particles """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import ManyAtomsSystem
-from mpliouville import MPLPDV
 import numpy as np
 from ase import Atoms
-from units import LJUnitsElement
 import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from rotagaporp.mpliouville import MPLPDV
+from utilities.units import LJUnitsElement
 
 argon_u = LJUnitsElement('Ar')
 symbols = 'ArAr'
@@ -18,18 +18,21 @@ atoms = Atoms(symbols, cell=cell, pbc=pbc, positions=positions, momenta=momenta)
 # numeric = {'tool': 'subs'}
 numeric = {'tool': 'theano'}
 # numeric = {'tool': 'lambdify', 'modules': ['math', 'numpy']}
+numeric = {'tool': 'lambdify', 'modules': 'tensorflow'}
 syst = ManyAtomsSystem(atoms, 'FFLennardJones', numeric=numeric)
 
 params = {'tstart': 0.0, 'tend': 14.0, 'tstep': 0.005, 'nterms': 4,
           'liouville': MPLPDV, 'tqdm': True, 'pbc': True}
-prop = Chebyshev(syst=syst, restart=True, restart_file='restart.pkl', **params)
-# prop = Chebyshev(syst=syst, **params)
-# prop.dump('restart.pkl')
-# exit()
+
+filename = 'restart.pkl'
+prop1 = Chebyshev(syst=syst, **params)
+prop1.dump(filename)
+
+prop2 = Chebyshev(syst=syst, restart=True, restart_file=filename, **params)
 print('propagation begins')
-prop.propagate()
-# prop.analyse()
-(ti, qt, pt) = prop.get_trajectory_3d()
+prop2.propagate()
+# prop2.analyse()
+(ti, qt, pt) = prop2.get_trajectory_3d()
 
 params = {'tstart': 0.0, 'tend': 14.0, 'tstep': 0.005, 'tqdm': True,
           'pbc': True}
diff --git a/classical/test_cheb_lj_2p_xd.py b/examples/test_cheb_lj_2p_xd.py
similarity index 86%
rename from classical/test_cheb_lj_2p_xd.py
rename to examples/test_cheb_lj_2p_xd.py
index 675879b9e2515bfa31e452364be395cb64f5ff3a..a39e07b7e8b3be6b6ac9fea3951533d17ea734e8 100644
--- a/classical/test_cheb_lj_2p_xd.py
+++ b/examples/test_cheb_lj_2p_xd.py
@@ -1,12 +1,12 @@
 """ Two particles in LJ potential - 1P, 1D and 3D implementation """
-from chebyshev import Chebyshev
-# from newtonian import Newtonian
-# from symplectic import Symplectic
-from systems import LennardJones, MPPS1D, MPPS3D
-from mpliouville import MPLPQP
-from liouville import LiouvillePowersQP
-import matplotlib.pyplot as plt
 import numpy as np
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+# from rotagaporp.newtonian import Newtonian
+# from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import LennardJones, MPPS1D, MPPS3D
+from rotagaporp.mpliouville import MPLPQP
+from rotagaporp.liouville import LiouvillePowersQP
 
 propagator = Chebyshev
 # propagator = Symplectic
@@ -19,11 +19,11 @@ params_1p = {'tstart': 0.0, 'tend': 2.0, 'tstep': 0.05, 'nterms': 4,
 
 M = [2.0, 2.0]
 syspar_3d = {'symbols': ['0x:z', '1x:z'], 'q0': [[0, 0, 0], [2, 0, 0]],
-          'p0': [[1, 0, 0], [-1, 0, 0]], 'masses': M}
+             'p0': [[1, 0, 0], [-1, 0, 0]], 'masses': M}
 ffpars_3d = [{'class': 'FFLennardJones', 'pair': ((0, 1, 2), (3, 4, 5))}]
 
 syspar_1d = {'symbols': ['0', '1'], 'q0': [[0], [2]], 'p0': [[1], [-1]],
-          'masses': M}
+             'masses': M}
 ffpars_1d = [{'class': 'FFLennardJones', 'pair': ((0,), (1,))}]
 
 syspar_1p = {
diff --git a/classical/test_cheb_lj_3p_3d.py b/examples/test_cheb_lj_3p_3d.py
similarity index 91%
rename from classical/test_cheb_lj_3p_3d.py
rename to examples/test_cheb_lj_3p_3d.py
index 72bdd74b31a52a6fc0c50a068e82a8498a860587..9652136d15176bf41bf24dc8fc9d84f76d4b28b5 100644
--- a/classical/test_cheb_lj_3p_3d.py
+++ b/examples/test_cheb_lj_3p_3d.py
@@ -1,13 +1,13 @@
 """ Three dimensional Lennard-Jones gas with three particles """
-from chebyshev import Chebyshev
-from symplectic import Symplectic
-from systems import ManyAtomsSystem
-from mpliouville import MPLPAQ
 import numpy as np
 from ase import Atoms
-from units import LJUnitsElement
-from trajectory import write_ase_trajectory
 import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from rotagaporp.mpliouville import MPLPAQ
+from utilities.units import LJUnitsElement
+from utilities.trajectory import write_ase_trajectory
 
 argon_u = LJUnitsElement('Ar')
 symbols = 'ArArAr'
diff --git a/classical/test_cheb_mo_2p_xd.py b/examples/test_cheb_mo_2p_xd.py
similarity index 89%
rename from classical/test_cheb_mo_2p_xd.py
rename to examples/test_cheb_mo_2p_xd.py
index 9ff27cebd60195bc3be6b352c402ff37a21dbea4..e4d32219f0490f4cf6ce682ba602a2736aaadfb6 100644
--- a/classical/test_cheb_mo_2p_xd.py
+++ b/examples/test_cheb_mo_2p_xd.py
@@ -1,12 +1,12 @@
 """ Two particles in Morse potential - 1P, 1D and 3D implementation """
-from chebyshev import Chebyshev
-# from newtonian import Newtonian
-# from symplectic import Symplectic
-from systems import Morse, MPPS1D, MPPS3D
-from mpliouville import MPLPQP
-from liouville import LiouvillePowersQP
-import matplotlib.pyplot as plt
 import numpy as np
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev
+# from rotagaporp.newtonian import Newtonian
+# from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import Morse, MPPS1D, MPPS3D
+from rotagaporp.mpliouville import MPLPQP
+from rotagaporp.liouville import LiouvillePowersQP
 
 propagator = Chebyshev
 # propagator = Symplectic
diff --git a/classical/test_cheb_recursion.py b/examples/test_cheb_recursion.py
similarity index 72%
rename from classical/test_cheb_recursion.py
rename to examples/test_cheb_recursion.py
index fae4d41b4a16f3633575d165b3d8614ccd4c73c5..006cf76473fca1088b89b1f6494440fc02274088 100644
--- a/classical/test_cheb_recursion.py
+++ b/examples/test_cheb_recursion.py
@@ -1,6 +1,8 @@
 """ Test Chebyshev recursion coefficients C^(n)_k """
-from combinatorics import cnk_chebyshev_T1, cnk_chebyshev_T2, cnk_chebyshev_T3
-from combinatorics import cnk_chebyshev_U1, cnk_chebyshev_U2, cnk_chebyshev_U3
+from rotagaporp.combinatorics import cnk_chebyshev_T1, cnk_chebyshev_T2
+from rotagaporp.combinatorics import cnk_chebyshev_T3
+from rotagaporp.combinatorics import cnk_chebyshev_U1, cnk_chebyshev_U2
+from rotagaporp.combinatorics import cnk_chebyshev_U3
 
 number = 32
 for signed in (True, False):
diff --git a/classical/test_ff_derivs.py b/examples/test_ff_derivs.py
similarity index 73%
rename from classical/test_ff_derivs.py
rename to examples/test_ff_derivs.py
index 6668b23c95b24aaa56000cd33b6b00884e419db4..24e2aa8316f2b3e31db65faae93054d40dfeb03b 100644
--- a/classical/test_ff_derivs.py
+++ b/examples/test_ff_derivs.py
@@ -1,19 +1,20 @@
+""" Test the force-field derivatives in 3D """
 import sympy as sp
 import numpy as np
-from systems import MPPS3D
+from rotagaporp.systems import MPPS3D
 
 syspar_3d = {'symbols': ['0x:z', '1x:z'], 'q0': [[1, 2, 3],
              [1.7071, 2.7071, 3]], 'p0': [[0.7071, 0.7071, 0], [0, 0, 0]],
              'masses':  np.full(2, 2.0)}
 
 ffpar_3d = {
-    'harm': [{'class': 'FFHarmonic', 'pair': ((0, 1, 2), (3, 4, 5)),
-             'params': {'distance': 3.5, 'kappa':  sp.Rational(1, 2)}}],
+    'harmonic': [{'class': 'FFHarmonic', 'pair': ((0, 1, 2), (3, 4, 5)),
+                  'params': {'distance': 3.5, 'kappa':  sp.Rational(1, 2)}}],
     'lj': [{'class': 'FFLennardJones', 'pair': ((0, 1, 2), (3, 4, 5)),
-           'params': {}}]
+            'params': {}}]
 }
 
-for ff in ('harm', 'lj'):
+for ff in ('harmonic', 'lj'):
     syst_3d = MPPS3D(ffpar_3d[ff], **syspar_3d)
     print(ff)
     print(list(syst_3d.ff[0].get_cart_derivs(2, explicit=False)))
diff --git a/classical/test_iln_00_1p_1d.py b/examples/test_iln_00_1p_1d.py
similarity index 71%
rename from classical/test_iln_00_1p_1d.py
rename to examples/test_iln_00_1p_1d.py
index d2336f78c8d4882a6ddee4038084d599992247d1..cab9bb134eaeb548a98b140be72d6218d6ab10df 100644
--- a/classical/test_iln_00_1p_1d.py
+++ b/examples/test_iln_00_1p_1d.py
@@ -1,4 +1,5 @@
-from liouville import q, p, LiouvillePowers
+""" Test the polynomial expansion for a free particle """
+from rotagaporp.liouville import q, p, LiouvillePowers
 
 order = 4
 ilnq = LiouvillePowers(order, function=q, hamiltonian=p**2/2,
diff --git a/classical/test_iln_1p.py b/examples/test_iln_1p.py
similarity index 95%
rename from classical/test_iln_1p.py
rename to examples/test_iln_1p.py
index 01cdb6dbde7e733b04c521fd9a295cf09db8b55e..f999bc433ad31e40c06df05af36a8a928a073e68 100644
--- a/classical/test_iln_1p.py
+++ b/examples/test_iln_1p.py
@@ -1,7 +1,6 @@
 """ Analysis of the terms of one-particle (iL)^n in one dimension """
-
 import sympy as sp
-from liouville import q, p, LiouvillePowers
+from rotagaporp.liouville import q, p, LiouvillePowers
 
 order = 5
 iln = LiouvillePowers(order, function=q, h_mixed=False, f_mixed=False)
diff --git a/classical/test_iln_2p_1d.py b/examples/test_iln_2p_1d.py
similarity index 95%
rename from classical/test_iln_2p_1d.py
rename to examples/test_iln_2p_1d.py
index a820b010f118b52a0dea0f8d6a2c5e3557aa3f08..4eb1e96c9802b7806e41ca91e4339411f4e6536b 100644
--- a/classical/test_iln_2p_1d.py
+++ b/examples/test_iln_2p_1d.py
@@ -1,7 +1,6 @@
 """ Analysis of the terms of two-particle (iL)^n in one dimension """
-
-from mpliouville import MPLP
 import sympy as sp
+from rotagaporp.mpliouville import MPLP
 
 nterms = 5
 Q = sp.Array(sp.symbols(['q1 q2'], real=True))
diff --git a/classical/test_iln_2p_1d_aq.py b/examples/test_iln_2p_1d_aq.py
similarity index 94%
rename from classical/test_iln_2p_1d_aq.py
rename to examples/test_iln_2p_1d_aq.py
index c43c11214c44777867d885a8dc0e7eeec773b0bd..6ac6535ab4c98081bc034d38f19d70104f249257 100644
--- a/classical/test_iln_2p_1d_aq.py
+++ b/examples/test_iln_2p_1d_aq.py
@@ -1,8 +1,7 @@
 """ Analysis of (iL)^n terms of a two-particle system in one dimension """
-
-from mpliouville import MPLPAQ, get_stats_expr, get_stats_repl
-import sympy as sp
 import yaml
+import sympy as sp
+from rotagaporp.mpliouville import MPLPAQ, get_stats_expr, get_stats_repl
 
 nterms = 4
 
diff --git a/classical/test_iln_2p_3d_dr.py b/examples/test_iln_2p_3d_dr.py
similarity index 90%
rename from classical/test_iln_2p_3d_dr.py
rename to examples/test_iln_2p_3d_dr.py
index 7825964fd53145cdab223030cb260689304a570e..53cbf787c0eccae51b05dacbb8fc7085aa15b491 100644
--- a/classical/test_iln_2p_3d_dr.py
+++ b/examples/test_iln_2p_3d_dr.py
@@ -1,9 +1,9 @@
 """ Analysis of the terms of three-particle (iL)^n in three dimensions """
-
-from mpliouville import MPLPDR, get_stats_expr, get_stats_repl
-from forcefield import FFList
-import sympy as sp
 import yaml
+import sympy as sp
+from rotagaporp.mpliouville import MPLPDR, get_stats_expr, get_stats_repl
+from rotagaporp.forcefield import FFList
+from rotagaporp import forcefield
 
 nterms = 3
 
diff --git a/classical/test_iln_2p_3d_dv.py b/examples/test_iln_2p_3d_dv.py
similarity index 89%
rename from classical/test_iln_2p_3d_dv.py
rename to examples/test_iln_2p_3d_dv.py
index ef3b3aa4c17528c075182fca4107fa3cd0e6db9d..dfad9c3cc64019ab5236b24f18346c555335951d 100644
--- a/classical/test_iln_2p_3d_dv.py
+++ b/examples/test_iln_2p_3d_dv.py
@@ -1,10 +1,9 @@
 """ Analysis of the terms of three-particle (iL)^n in three dimensions """
-
-from mpliouville import MPLPDV, get_stats_expr
-from systems import MPPS3D
+import yaml
 import sympy as sp
 import numpy as np
-import yaml
+from rotagaporp.mpliouville import MPLPDV, get_stats_expr
+from rotagaporp.systems import MPPS3D
 
 nterms = 3
 
@@ -21,7 +20,7 @@ system = MPPS3D(ffparams, symbols, masses, qval, pval, pbc=pbc, cell=cell,
                 numeric=numeric)
 
 iln = MPLPDV(nterms, system, variant='r')
-print(iln.get_val_float(qval, pval))
+print(iln.get_val_float(system.q0, system.p0))
 print(yaml.dump(get_stats_expr(iln.ilnqexpr, details=True)))
 print(yaml.dump(get_stats_expr(iln.ilnpexpr, details=True)))
 
diff --git a/classical/test_iln_2p_3d_nr_cheb.py b/examples/test_iln_2p_3d_nr_cheb.py
similarity index 79%
rename from classical/test_iln_2p_3d_nr_cheb.py
rename to examples/test_iln_2p_3d_nr_cheb.py
index babc21d5042759c62a42fe002dba36859772a12b..e09e6bd0abce200340941ff91261302617f5cf9a 100644
--- a/classical/test_iln_2p_3d_nr_cheb.py
+++ b/examples/test_iln_2p_3d_nr_cheb.py
@@ -1,17 +1,16 @@
 """ Non-recursive Chebyshev expansion of (iL)^n for two particles in 3D """
-
-from mpliouville import MPLPDV, MPLPNR
-from systems import MPPS3D
-import numpy as np
 import time
+import numpy as np
+from rotagaporp.mpliouville import MPLPDV, MPLPNR
+from rotagaporp.systems import MPPS3D
 
 nterms = 3
 
 ffparams = [{'class': 'FFLennardJones', 'pair': ((0, 1, 2), (3, 4, 5))}]
 symbols = ['0x:z', '1x:z']
-masses = np.array([1, 1, 1, 1, 1, 1])
-qval = np.random.rand(6)
-pval = np.random.rand(6)
+masses = np.array([1, 1])
+qval = np.random.rand(2, 3)
+pval = np.random.rand(2, 3)
 # numeric={'tool': 'theano'}
 # numeric={'tool': 'lambdify', 'modules': 'math'}
 numeric={'tool': 'subs'}
@@ -32,8 +31,8 @@ ilnnr = MPLPNR(nterms, system, recursion='chebyshev', sign=1)
 te = time.time()
 print("preparation time for MPLPNR:", te-ts)
 
-assert np.allclose(ilnnr.get_val_float(qval, pval),
-                   ilndv.get_val_float(qval, pval))
+assert np.allclose(ilnnr.get_val_float(system.q0, system.p0),
+                   ilndv.get_val_float(system.q0, system.p0))
 
 qarr = np.random.rand(1000, 6)
 parr = np.random.rand(1000, 6)
diff --git a/classical/test_iln_2p_3d_nr_newt.py b/examples/test_iln_2p_3d_nr_newt.py
similarity index 79%
rename from classical/test_iln_2p_3d_nr_newt.py
rename to examples/test_iln_2p_3d_nr_newt.py
index 38415a5f59154bb5a3387b7f748b9ee8cdd99790..bacd98346727a10fa63c1f1381cc3d04bd33f51f 100644
--- a/classical/test_iln_2p_3d_nr_newt.py
+++ b/examples/test_iln_2p_3d_nr_newt.py
@@ -1,18 +1,17 @@
 """ Non-recursive Newtonian expansion of (iL)^n for two particles in 3D """
-
-from mpliouville import MPLPDV, MPLPNR
-from newtonian import leja_points_ga
-from systems import MPPS3D
-import numpy as np
 import time
+import numpy as np
+from rotagaporp.mpliouville import MPLPDV, MPLPNR
+from rotagaporp.newtonian import leja_points_ga
+from rotagaporp.systems import MPPS3D
 
 nterms = 3
 
 ffparams = [{'class': 'FFLennardJones', 'pair': ((0, 1, 2), (3, 4, 5))}]
 symbols = ['0x:z', '1x:z']
-masses = np.array([1, 1, 1, 1, 1, 1])
-qval = np.random.rand(6)
-pval = np.random.rand(6)
+masses = np.array([1, 1])
+qval = np.random.rand(2, 3)
+pval = np.random.rand(2, 3)
 numeric={'tool': 'theano'}
 # numeric={'tool': 'lambdify', 'modules': 'math'}
 # numeric={'tool': 'subs'}
@@ -38,8 +37,8 @@ ilnnr = MPLPNR(nterms, system, recursion='newton', nodes=lambdas,
 te = time.time()
 print("preparation time for MPLPNR:", te-ts)
 
-assert np.allclose(ilnnr.get_val_float(qval, pval),
-                   ilndv.get_val_float(qval, pval))
+assert np.allclose(ilnnr.get_val_float(system.q0, system.p0),
+                   ilndv.get_val_float(system.q0, system.p0))
 
 qarr = np.random.rand(1000, 6)
 parr = np.random.rand(1000, 6)
diff --git a/classical/test_iln_3p_1d_aq.py b/examples/test_iln_3p_1d_aq.py
similarity index 92%
rename from classical/test_iln_3p_1d_aq.py
rename to examples/test_iln_3p_1d_aq.py
index f41e71df47aa435358a64ab6c11463e3736b2c9c..49e18ea29f9f1a1427a3fc6e586203a8d8162699 100644
--- a/classical/test_iln_3p_1d_aq.py
+++ b/examples/test_iln_3p_1d_aq.py
@@ -1,9 +1,8 @@
 """ Analysis of the terms of three-particle (iL)^n in one dimension """
-
-from mpliouville import MPLPAQ, get_stats_expr, get_stats_repl
-from forcefield import FFList
-import sympy as sp
 import yaml
+import sympy as sp
+from rotagaporp.mpliouville import MPLPAQ, get_stats_expr, get_stats_repl
+from rotagaporp.forcefield import FFList
 
 nterms = 3
 
diff --git a/classical/test_iln_3p_1d_vq.py b/examples/test_iln_3p_1d_vq.py
similarity index 92%
rename from classical/test_iln_3p_1d_vq.py
rename to examples/test_iln_3p_1d_vq.py
index def7e051bbdf932c6b714ae1f14746edcbef1e54..f76ff974ac534038de89b9f392413b003f43272d 100644
--- a/classical/test_iln_3p_1d_vq.py
+++ b/examples/test_iln_3p_1d_vq.py
@@ -1,9 +1,8 @@
 """ Analysis of the terms of three-particle (iL)^n in one dimension """
-
-from mpliouville import MPLPVQ, get_stats_expr, get_stats_repl
-from forcefield import FFList
-import sympy as sp
 import yaml
+import sympy as sp
+from rotagaporp.mpliouville import MPLPVQ, get_stats_expr, get_stats_repl
+from rotagaporp.forcefield import FFList
 
 nterms = 3
 
diff --git a/classical/test_iln_4p_3d.py b/examples/test_iln_4p_3d.py
similarity index 77%
rename from classical/test_iln_4p_3d.py
rename to examples/test_iln_4p_3d.py
index ec54a8e0f30fb02918b738daf754c8afbd725d82..bcaea254fe05b54a2f2d48200395e78bcdc49dd1 100644
--- a/classical/test_iln_4p_3d.py
+++ b/examples/test_iln_4p_3d.py
@@ -1,11 +1,10 @@
 """ Analysis of the terms of four-particle (iL)^n in three dimensions """
-
-from itertools import combinations
-import numpy as np
-from mpliouville import MPLPDV, get_stats_expr, get_stats_repl
-from systems import MPPS3D
 import yaml
 import time
+import numpy as np
+from itertools import combinations
+from rotagaporp.mpliouville import MPLPDV, get_stats_expr, get_stats_repl
+from rotagaporp.systems import MPPS3D
 
 nprtcls = 5
 numeric = {'tool': 'subs'}
@@ -14,9 +13,9 @@ symbols = [str(i)+'x:z' for i in range(nprtcls)]
 indic = list(range(3*nprtcls))
 pairs = list(combinations(zip(indic[0::3], indic[1::3], indic[2::3]), 2))
 ffparams = [{'class': 'FFLennardJones', 'pair': p, 'params': {}} for p in pairs]
-masses = np.full((nprtcls*3, ), 1)
-qval = np.full((nprtcls*3, ), 1) # these are not valid
-pval = np.full((nprtcls*3, ), 0)
+masses = np.full((nprtcls, ), 1)
+qval = np.full((nprtcls, 3), 1) # these are not valid
+pval = np.full((nprtcls, 3), 0)
 system = MPPS3D(ffparams, symbols, masses, qval, pval, numeric=numeric)
 
 nterms = 3
diff --git a/examples/test_iln_np_3d.py b/examples/test_iln_np_3d.py
new file mode 100644
index 0000000000000000000000000000000000000000..0c116ddb1723ef152aabb58804d8e9d13d625966
--- /dev/null
+++ b/examples/test_iln_np_3d.py
@@ -0,0 +1,27 @@
+""" Analysis of the terms of three-particle (iL)^n in three dimensions """
+import yaml
+import numpy as np
+from itertools import combinations
+from rotagaporp.mpliouville import MPLPDV, get_stats_expr, get_stats_repl
+from rotagaporp.systems import MPPS3D
+
+nprtcls = 3
+# numeric = {'tool': 'subs'}
+numeric = {'tool': 'lambdify', 'modules': 'math'}
+
+symbols = [str(i)+'x:z' for i in range(nprtcls)]
+indic = list(range(3*nprtcls))
+pairs = list(combinations(zip(indic[0::3], indic[1::3], indic[2::3]), 2))
+ffparams = [{'class': 'FFLennardJones', 'pair': p, 'params': {}} for p in pairs]
+masses = np.full((nprtcls, ), 1)
+qval = np.random.rand(nprtcls, 3)
+pval = np.random.rand(nprtcls, 3)
+system = MPPS3D(ffparams, symbols, masses, qval, pval, numeric=numeric)
+
+nterms = 4
+
+iln = MPLPDV(nterms, system, variant='r', numeric=numeric)
+# print(yaml.dump(get_stats_expr(iln.ilnqexpr, details=True)))
+# print(yaml.dump(get_stats_expr(iln.ilnpexpr, details=True)))
+# print(yaml.dump(get_stats_repl(iln.repl, details=True)))
+print(iln.get_val_float(system.q0, system.p0))
diff --git a/examples/test_liouville_qp_ev_1p.py b/examples/test_liouville_qp_ev_1p.py
new file mode 100644
index 0000000000000000000000000000000000000000..3d0a89a64a4140c59ea64f61f5f419ece5ad1af6
--- /dev/null
+++ b/examples/test_liouville_qp_ev_1p.py
@@ -0,0 +1,59 @@
+""" Numerical determination of Liouville eigenvalues for one particle """
+import numpy as np
+import matplotlib
+matplotlib.use('Qt5Cairo')
+import matplotlib.pyplot as plt
+from rotagaporp.systems import Morse, Harmonic, Anharmonic, LennardJones
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.chebyshev import Chebyshev
+
+params = {'tstart': 0., 'tend': 100., 'tstep': 0.01, 'delta': 1., 'nterms': 5}
+syspar = {'q0': 3., 'p0': 0.}
+syst = Morse(**syspar)
+# syst = Harmonic(**syspar)
+# syst = Anharmonic(**syspar)
+# syst = LennardJones(**syspar)
+prop = Chebyshev(syst=syst, **params)
+# prop = Symplectic(syst=syst, **params)
+prop.propagate()
+# prop.analyse()
+(ti, qt, pt) = prop.get_trajectory()
+
+half = len(ti)//2
+
+fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
+ax1.plot(ti, qt, label='q')
+ax1.plot(ti, pt, label='p')
+ax1.set_xlabel(r'$t$')
+ax1.set_ylabel(r'$q,\quad p$')
+ax1.legend()
+
+ax2.plot(qt, pt)
+ax2.set_xlabel(r'$q$')
+ax2.set_ylabel(r'$p$')
+
+variable = pt
+# '''
+sigma = (params['tend']-params['tstart'])/3
+center = (params['tstart']+params['tend'])/2
+corrf = np.correlate(variable, variable, mode='same')
+gauss = np.exp(-(ti-center)**2/(2*sigma**2))/(np.sqrt(2*np.pi)*sigma)
+ax3.plot(ti, corrf*gauss)
+ax3.set_xlabel(r'$\tau$')
+ax3.set_ylabel(r'$C(\tau)=\int_{-\infty}^{+\infty}\ p(t+\tau)p(t)dt$')
+fftf = np.fft.fft(corrf*gauss)
+'''
+cfft = np.fft.fft(variable)
+fftf = cfft*np.conjugate(cfft)
+'''
+freq = np.fft.fftfreq(len(ti), d=params['tstep'])*2*np.pi
+
+# ax4.plot(freq, fftf.real, label='Re(F)')
+# ax4.plot(freq, fftf.imag, label='Im(F)')
+ax4.plot(freq[:half], np.absolute(fftf[:half]), label=r'$|F(\omega)|$')
+ax4.set_xlabel(r'$\omega$')
+ax4.set_ylabel(r'$F(\omega)={\cal F} C(\tau)$')
+ax4.legend()
+
+fig.subplots_adjust(hspace=0.3)
+plt.show()
diff --git a/examples/test_liouville_rho_ev_1p.py b/examples/test_liouville_rho_ev_1p.py
new file mode 100644
index 0000000000000000000000000000000000000000..e12ba59d28f5a2c445fd44a59cc45f90cb456b7d
--- /dev/null
+++ b/examples/test_liouville_rho_ev_1p.py
@@ -0,0 +1,79 @@
+""" Numerical determination of Liouville eigenfunctions for one particle """
+import numpy as np
+import matplotlib
+matplotlib.use('Qt5Cairo')
+import matplotlib.pyplot as plt
+from rotagaporp.systems import Morse, Harmonic, Anharmonic, LennardJones
+from rotagaporp.symplectic import Symplectic
+
+params = {'tstart': 0., 'tend': 1000., 'tstep': 0.1}
+syspar = {'q0': 2., 'p0': 0.0}
+
+syst = Morse(**syspar)
+# syst = Harmonic(**syspar)
+# syst = Anharmonic(**syspar)
+# syst = LennardJones(**syspar)
+prop = Symplectic(syst=syst, **params)
+prop.propagate()
+(ti, qt, pt) = prop.get_trajectory()
+half = len(ti)//2
+
+fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
+ax1.plot(ti, qt, label='q')
+ax1.plot(ti, pt, label='p')
+ax1.set_xlabel(r'$t$')
+ax1.set_ylabel(r'$q,\quad p$')
+ax1.legend()
+
+ax2.plot(qt, pt)
+ax2.set_xlabel(r'$q$')
+ax2.set_ylabel(r'$p$')
+
+border = 0.05 # 5%
+pmin = np.min(pt)*(1-border)
+pmax = np.max(pt)*(1+border)
+qmin = np.min(qt)*(1-border)
+qmax = np.max(qt)*(1+border)
+qsigma = (qmax-qmin)/10
+psigma = (pmax-pmin)/10
+prange = np.linspace(pmin, pmax, 80)
+qrange = np.linspace(qmin, qmax, 80)
+qmesh, pmesh = np.meshgrid(qrange, prange, indexing='ij')
+rho = np.empty((len(ti), len(qrange), len(prange)), dtype=np.float)
+for tind in range(len(ti)):
+    rho[tind, :, :] = (np.exp(-(qt[tind]-qmesh)**2/(2*qsigma**2)
+                              -(pt[tind]-pmesh)**2/(2*psigma**2))
+                        /(2*np.pi)*qsigma*psigma)
+
+# '''
+cfft = np.fft.fft(rho, axis=0)
+fftf = cfft*np.conjugate(cfft)
+'''
+corrf = np.empty(rho.shape, dtype=np.float)
+sigma = (params['tend']-params['tstart'])/3
+center = (params['tstart']+params['tend'])/2
+for i in range(rho.shape[1]):
+    for j in range(rho.shape[2]):
+        corrf[:, i, j] = np.correlate(rho[:, i, j], rho[:, i, j], mode='same')
+fftf = np.empty(rho.shape, dtype=np.complex)
+for i in range(rho.shape[1]):
+    for j in range(rho.shape[2]):
+        fftf[:, i, j] = np.fft.fft(corrf[:, i, j])
+'''
+freq = np.fft.fftfreq(len(ti), d=params['tstep'])*2*np.pi
+
+fftcjg = fftf*np.conjugate(fftf)
+avfftf = np.trapz(np.trapz(fftcjg, x=prange, axis=2), x=qrange, axis=1)
+ax3.plot(freq[:half], avfftf[:half].real, label='integrated')
+ax3.set_xlabel(r'$\omega$')
+ax3.set_ylabel(r'$F(\omega)={\cal F} C(\tau)$')
+ax3.legend()
+
+show_ev = 1.4
+show_ev_index = np.min(np.where(freq >= show_ev))
+ax4.contourf(qmesh, pmesh, fftf[show_ev_index, :, :].real)
+ax4.set_xlabel('q')
+ax4.set_ylabel('p')
+
+fig.subplots_adjust(hspace=0.3)
+plt.show()
diff --git a/examples/test_verlet_00_2p_3d_pbc.py b/examples/test_verlet_00_2p_3d_pbc.py
new file mode 100644
index 0000000000000000000000000000000000000000..0f7890a3cee76f7b487feb373c8ec1bc737e7d34
--- /dev/null
+++ b/examples/test_verlet_00_2p_3d_pbc.py
@@ -0,0 +1,59 @@
+""" Three dimensional Lennard-Jones gas with two particles """
+import numpy as np
+from ase import Atoms
+import matplotlib.pyplot as plt
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from utilities.units import LJUnitsElement
+from utilities.trajectory import write_ase_trajectory
+
+argon_u = LJUnitsElement('Ar')
+symbols = 'ArAr'
+pbc = (True, True, True)
+cell = np.array([15, 15, 15])*argon_u.sigmaA
+positions = np.array([[1, 0, 0], [3, 0, 0]])*argon_u.sigmaA
+momenta = np.array([[0, 0, 0], [-1, 0, 0]])*argon_u.momentumA
+atoms = Atoms(symbols, cell=cell, pbc=pbc, positions=positions, momenta=momenta)
+numeric = {'tool': 'lambdify', 'modules': ['math', 'numpy']}
+# numeric = {'tool': 'theano'}
+# numeric = {'tool': 'subs'}
+syst = ManyAtomsSystem(atoms, 'FFZero', numeric=numeric)
+
+params = {'tstart': 0.0, 'tend': 28.0, 'tstep': 0.005, 'tqdm': True, 'pbc': True}
+prop_ref = Symplectic(syst=syst, **params)
+prop_ref.propagate()
+prop_ref.analyse()
+enr, err = prop_ref.en, prop_ref.er
+(tir, qtr, ptr) = prop_ref.get_trajectory_3d(step=1)
+write_ase_trajectory(prop_ref, atoms, step=1)
+
+# distance between the two particles with time
+rtr = [syst.get_distances(q)[0] for q in qtr]
+
+fig = plt.figure()
+plot = fig.add_subplot(411)
+plt.plot(tir, qtr[:, 0, 0], label='q0 verlet')
+plt.plot(tir, qtr[:, 1, 0], label='q1 verlet')
+plot.set_ylabel(r'$q_0,\quad q_1$')
+plt.legend()
+
+plot = fig.add_subplot(412)
+plt.plot(tir, ptr[:, 0, 0], label='p0 verlet')
+plt.plot(tir, ptr[:, 1, 0], label='p1 verlet')
+plot.set_ylabel(r'$p_0,\quad p_1$')
+plt.legend()
+
+plot = fig.add_subplot(413)
+plt.plot(tir, rtr, label='r verlet')
+plot.set_ylabel(r'$r(t)$')
+plot.set_xlabel('time')
+plt.legend()
+
+plot = fig.add_subplot(414)
+plt.plot(tir, enr, label='energy, verlet')
+plt.plot(tir, err, label='energy error, verlet')
+plot.set_ylabel(r'$E(r)\quad\Delta E(t)$')
+plot.set_xlabel('time')
+plt.legend()
+
+plt.show()
diff --git a/examples/test_verlet_argon.py b/examples/test_verlet_argon.py
new file mode 100644
index 0000000000000000000000000000000000000000..413ffb92c89ff93393c8bd5d1a43be0e5d976e91
--- /dev/null
+++ b/examples/test_verlet_argon.py
@@ -0,0 +1,39 @@
+""" Molecular dynamics of argon at constant energy """
+import time
+from ase import units
+from ase.build import bulk
+from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
+from rotagaporp.systems import ManyAtomsSystem
+from rotagaporp.symplectic import Symplectic
+from utilities.trajectory import write_ase_trajectory
+
+T0_lj = 3.0 # initial temperature in Ar LJ units
+rho_lj = 0.1 # 0.88 # density in Ar LJ units
+sigma_a = 3.405*units.Angstrom # in Angstrom for Ar
+epsilon_a = 1.65e-21*units.J # in Joule for Ar
+T = T0_lj * epsilon_a / units.kB # initial temperature in Kelvin
+rho = rho_lj/sigma_a**3 # density in atoms per cubic Angstrom
+n_cell = 2 # number of unit cells in one direction
+n_atoms = 4*n_cell**3 # total number of atoms, 4 atoms per unit cell
+l_box = (n_atoms/rho)**(1./3) # edge length of simulation box
+l_cell = l_box / n_cell # latice parameter
+
+atoms = bulk('Ar', 'fcc', a=l_cell, cubic=True).repeat((n_cell, n_cell, n_cell))
+rho_real = sum(atoms.get_masses())/atoms.get_volume()*units.m**3/units.kg
+print('density, kg/m^3', rho_real)
+
+prep_start = time.time()
+MaxwellBoltzmannDistribution(atoms, T*units.kB)
+syst = ManyAtomsSystem(atoms, 'FFLennardJones', numeric={'tool': 'theano'})
+params = {'tstart': 0.0, 'tend': 10.0, 'tstep': 0.001, 'tqdm': True, 'pbc': True}
+prop = Symplectic(syst=syst, **params)
+prep_end = time.time()
+
+print('preparation time:', prep_end-prep_start)
+print('propagation begins')
+prop_start = time.time()
+prop.propagate()
+prop_end = time.time()
+print('propagation time:', prop_end-prop_start)
+prop.analyse(step=10)
+write_ase_trajectory(prop, atoms, step=10)
diff --git a/classical/test_verlet_lj_2p_3d.py b/examples/test_verlet_lj_2p_3d.py
similarity index 89%
rename from classical/test_verlet_lj_2p_3d.py
rename to examples/test_verlet_lj_2p_3d.py
index 55162df08e3d38a741a2fc049037ed3c91ebcc37..c830d5d73002fec0618cfd5e2f3fb4376f4123ea 100644
--- a/classical/test_verlet_lj_2p_3d.py
+++ b/examples/test_verlet_lj_2p_3d.py
@@ -1,11 +1,11 @@
 """ Three dimensional Lennard-Jones gas with two particles """
-from symplectic import Symplectic
-from systems import ManyAtomsSystem
 import numpy as np
 from ase import Atoms
-from units import LJUnitsElement
-from trajectory import write_ase_trajectory
 import matplotlib.pyplot as plt
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from utilities.units import LJUnitsElement
+from utilities.trajectory import write_ase_trajectory
 
 argon_u = LJUnitsElement('Ar')
 symbols = 'ArAr'
diff --git a/classical/test_verlet_lj_2p_3d_pbc.py b/examples/test_verlet_lj_2p_3d_pbc.py
similarity index 90%
rename from classical/test_verlet_lj_2p_3d_pbc.py
rename to examples/test_verlet_lj_2p_3d_pbc.py
index a68be018a7870ef68f57ae5fb027590ba1acf361..8319e90f160b92b1b526516aa9f7116be02579f9 100644
--- a/classical/test_verlet_lj_2p_3d_pbc.py
+++ b/examples/test_verlet_lj_2p_3d_pbc.py
@@ -1,11 +1,11 @@
 """ Three dimensional Lennard-Jones gas with two particles """
-from symplectic import Symplectic
-from systems import ManyAtomsSystem
 import numpy as np
 from ase import Atoms
-from units import LJUnitsElement
-from trajectory import write_ase_trajectory
 import matplotlib.pyplot as plt
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from utilities.units import LJUnitsElement
+from utilities.trajectory import write_ase_trajectory
 
 argon_u = LJUnitsElement('Ar')
 symbols = 'ArAr'
diff --git a/classical/test_verlet_lj_3p_3d_pbc.py b/examples/test_verlet_lj_3p_3d_pbc.py
similarity index 92%
rename from classical/test_verlet_lj_3p_3d_pbc.py
rename to examples/test_verlet_lj_3p_3d_pbc.py
index c1750241386b004177ff91946a83fbee42a56c74..20b227e29d0b12f7aca3cdfa88de5ff63e05249a 100644
--- a/classical/test_verlet_lj_3p_3d_pbc.py
+++ b/examples/test_verlet_lj_3p_3d_pbc.py
@@ -1,11 +1,11 @@
 """ Three dimensional Lennard-Jones gas with three particles """
-from symplectic import Symplectic
-from systems import ManyAtomsSystem
 import numpy as np
 from ase import Atoms
-from units import LJUnitsElement
-from trajectory import write_ase_trajectory
 import matplotlib.pyplot as plt
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import ManyAtomsSystem
+from utilities.units import LJUnitsElement
+from utilities.trajectory import write_ase_trajectory
 
 argon_u = LJUnitsElement('Ar')
 symbols = 'ArArAr'
diff --git a/classical/time_1p_nr_cheb.py b/examples/time_1p_nr_cheb.py
similarity index 87%
rename from classical/time_1p_nr_cheb.py
rename to examples/time_1p_nr_cheb.py
index e0f4c03922200eebf39f397de106780bbcb988ee..6d3d24e47d04cfaa0756fd1a62d0321411d54a8a 100644
--- a/classical/time_1p_nr_cheb.py
+++ b/examples/time_1p_nr_cheb.py
@@ -1,10 +1,9 @@
 """ Test and time non-recursive Chebyshev expansion for one particle in 1D """
-
+import time
 import sympy as sp
 import numpy as np
-from liouville import LiouvillePowersNR, LiouvillePowersQP
-from systems import Harmonic, Anharmonic, Morse, LennardJones
-import time
+from rotagaporp.liouville import LiouvillePowersNR, LiouvillePowersQP
+from rotagaporp.systems import Harmonic, Anharmonic, Morse, LennardJones
 
 nterms = 16
 syst = LennardJones()
diff --git a/classical/time_1p_nr_newt.py b/examples/time_1p_nr_newt.py
similarity index 86%
rename from classical/time_1p_nr_newt.py
rename to examples/time_1p_nr_newt.py
index 0fe707a5df96370866cf057c3c6d8dae59a2a29c..52db0b8db37eea7da96898ef0a38979512f20373 100644
--- a/classical/time_1p_nr_newt.py
+++ b/examples/time_1p_nr_newt.py
@@ -1,11 +1,10 @@
 """ Test non-recursive Newton interpolation for one particle in 1D """
-
+import time
 import sympy as sp
 import numpy as np
-from liouville import LiouvillePowersNR, LiouvillePowersQP
-from systems import Harmonic, Anharmonic, Morse, LennardJones
-from newtonian import leja_points_ga
-import time
+from rotagaporp.liouville import LiouvillePowersNR, LiouvillePowersQP
+from rotagaporp.systems import Harmonic, Anharmonic, Morse, LennardJones
+from rotagaporp.newtonian import leja_points_ga
 
 nterms = 3
 syst = LennardJones()
diff --git a/examples/time_cheb_lj_np_3d_theano.py b/examples/time_cheb_lj_np_3d_theano.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac210c866c18383d3cc2659e53c3901d33b9ba4d
--- /dev/null
+++ b/examples/time_cheb_lj_np_3d_theano.py
@@ -0,0 +1,77 @@
+""" Scalability measurement of the construction of (iL)^n f expressions for an
+    N-particle 3D system """
+import time
+import traceback
+import numpy as np
+from itertools import combinations
+from pandas.io.json import json_normalize
+from rotagaporp.mpliouville import MPLPAQ, MPLPVQ, MPLPDR, MPLPDV, MPLPNR
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.systems import MPPS3D
+
+cases = (
+     (MPLPAQ, None, {'tool': 'theano'}),
+#     (MPLPVQ, None, {'tool': 'theano'}),
+#     (MPLPDR, 'r', {'tool': 'theano'}),
+#     (MPLPDR, 'dr', {'tool': 'theano'}),
+     (MPLPDV, 'r', {'tool': 'theano'}),
+     (MPLPDV, 'dv', {'tool': 'theano'}),
+     (MPLPNR, 'r', {'tool': 'theano'}),
+     (MPLPNR, 'dv', {'tool': 'theano'}),
+)
+
+params = {'tstart': 0.0, 'tend': 10.0, 'tstep': 0.001}
+qvals = np.array([[0, 0, 0], [2, 0, 0], [0, 2, 0],
+                  [0, 0, 2], [2, 2, 0], [0, 2, 2]])
+pvals = np.array([[0, 0, 0], [-1, 0, 0], [0, -1, 0],
+                  [0, 0, -1], [-1, -1, 0], [0, -1, -1]])
+
+statistics = []
+for nprtcls in range(2, 6):
+    print('nprtcls:', nprtcls)
+    symbols = [str(i)+'x:z' for i in range(nprtcls)]
+    indic = list(range(3*nprtcls))
+    pairs = list(combinations(zip(indic[0::3], indic[1::3], indic[2::3]), 2))
+    ffparams = [{'class': 'FFLennardJones', 'pair': p, 'params': {}} for p in pairs]
+    masses = np.full((nprtcls, ), 1)
+    qval = qvals[:nprtcls]
+    pval = pvals[:nprtcls]
+
+    for nterms in range(2, 6):
+        print('nterms:', nterms)
+        for mplp, variant, numeric in cases:
+            print('class:', mplp.__name__, 'variant:', variant,
+                  'numeric:', numeric)
+            stats = {}
+            stats['propagator'] = 'chebyshev'
+            stats['class'] = mplp.__name__
+            stats['variant'] = variant
+            stats['nterms'] = nterms
+            stats['numeric tool'] = numeric['tool']
+            stats['nparticles'] = nprtcls
+            cheb_params = params.copy()
+            cheb_params['nterms'] = nterms
+            cheb_params['liouville'] = mplp
+            cheb_params['syst'] = MPPS3D(ffparams, symbols, masses, qval, pval,
+                                         numeric=numeric)
+            cheb_params['liou_params'] = {'numeric': numeric}
+            if variant is not None:
+                cheb_params['liou_params']['variant'] = variant
+            try:
+                ts = time.time()
+                cheb = Chebyshev(**cheb_params)
+                te = time.time()
+                stats['preparation time'] = te - ts
+            except Exception as err:
+                traceback.print_exc()
+            try:
+                ts = time.time()
+                cheb.propagate()
+                te = time.time()
+                stats['propagation time'] = te - ts
+            except Exception as err:
+                traceback.print_exc()
+            statistics.append(stats)
+
+df = json_normalize(statistics)
+df.to_json('time_cheb_lj_np_3d_theano.json', double_precision=15)
diff --git a/classical/time_pbc.py b/examples/time_pbc.py
similarity index 81%
rename from classical/time_pbc.py
rename to examples/time_pbc.py
index 4e21c2dbc687a7b2d1366e18fb62f46947e36637..3369c762e12c4e3b2d90a751b68138211c9ba3ef 100644
--- a/classical/time_pbc.py
+++ b/examples/time_pbc.py
@@ -1,10 +1,11 @@
+""" Performance of different PBC algorithms """
 import numpy as np
-from timing import timeit
-from pbc import get_distance_matrix_lc
-from pbc import get_distance_matrix_bc
-from pbc import wrap_positions as wrap1
 from ase.geometry import get_distances
 from ase.geometry import wrap_positions as wrap2
+from rotagaporp.pbc import get_distance_matrix_lc
+from rotagaporp.pbc import get_distance_matrix_bc
+from rotagaporp.pbc import wrap_positions as wrap1
+from utilities.timing import timeit
 
 @timeit
 def time_pbc1(posvec, cell):
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ade80f932df56c6e0d14d7723f2f742a5a46a2de
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,10 @@
+ase==3.17.0
+cloudpickle==0.8.1
+matplotlib==3.0.3
+mpmath==1.1.0
+nose==1.3.7
+numpy==1.16.2
+scipy==1.2.1
+sympy==1.3
+Theano==1.0.4
+tqdm==4.31.1
diff --git a/rotagaporp/.gitignore b/rotagaporp/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..bee8a64b79a99590d5303307144172cfe824fbf7
--- /dev/null
+++ b/rotagaporp/.gitignore
@@ -0,0 +1 @@
+__pycache__
diff --git a/rotagaporp/__init__.py b/rotagaporp/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..659f3d97f3fb09fd455bc6d26a81f38eadbd0fa1
--- /dev/null
+++ b/rotagaporp/__init__.py
@@ -0,0 +1,2 @@
+import os, sys
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
diff --git a/classical/anharmonic.py b/rotagaporp/anharmonic.py
similarity index 100%
rename from classical/anharmonic.py
rename to rotagaporp/anharmonic.py
diff --git a/classical/ballistic.py b/rotagaporp/ballistic.py
similarity index 100%
rename from classical/ballistic.py
rename to rotagaporp/ballistic.py
diff --git a/classical/chebyshev.py b/rotagaporp/chebyshev.py
similarity index 100%
rename from classical/chebyshev.py
rename to rotagaporp/chebyshev.py
diff --git a/classical/combinatorics.py b/rotagaporp/combinatorics.py
similarity index 100%
rename from classical/combinatorics.py
rename to rotagaporp/combinatorics.py
diff --git a/classical/forcefield.py b/rotagaporp/forcefield.py
similarity index 98%
rename from classical/forcefield.py
rename to rotagaporp/forcefield.py
index 075374452a2042c2a54222c4a9e4875da325f043..c1751fbd8eaa63732b323eb62775ec5427fa0f9a 100644
--- a/classical/forcefield.py
+++ b/rotagaporp/forcefield.py
@@ -134,6 +134,8 @@ class FFList(list):
     def __init__(self, ffparams, Q):
         """ ffparams: [{'class': str, 'pair': (), 'params': {}}] """
         namespace = __import__(__name__)
+        for subpackage in __name__.split('.')[1:]:
+            namespace = getattr(namespace, subpackage)
         for index, ffpar in enumerate(ffparams):
             params = ffpar.get('params', {})
             params['label'] = str(index)
diff --git a/classical/generic.py b/rotagaporp/generic.py
similarity index 82%
rename from classical/generic.py
rename to rotagaporp/generic.py
index f949ae01c89dc8808f995ceee0adc48fb17c5d47..aec0e73a9f5347072ccb535894ece7607ac6edf3 100644
--- a/classical/generic.py
+++ b/rotagaporp/generic.py
@@ -12,39 +12,34 @@ def _float(fpprec, f):
 
 def spfloat(f, prec):
     """ creates a sympy float with specified precision from a python float """
-    from decimal import Decimal
     if isinstance(f, float):
-        x, y = Decimal(repr(f)).as_integer_ratio()
-        return sp.N(x, prec)/sp.N(y, prec)
-    else:
-        return f
+        from decimal import Decimal
+        from fractions import Fraction
+        # only in Python 3.6 or later
+        # x, y = Decimal(repr(f)).as_integer_ratio()
+        # ret = sp.N(x, prec)/sp.N(y, prec)
+        frac = Fraction.from_decimal(Decimal(repr(f)))
+        ret = sp.N(frac.numerator, prec)/sp.N(frac.denominator, prec)
+    return ret
 
 def flatten(x):
     """ Return a flat list from a deeply nested iterables """
     from collections import Iterable
     return (a for i in x for a in flatten(i)) if isinstance(x, Iterable) else [x]
 
-def uniq1(seq):
-    """ filter out duplicated elements and return the unique ones, faster """
+def uniq(seq):
+    """ filter out duplicated elements and return the unique ones """
     seen = set()
     seen_add = seen.add
     return [x for x in seq if not (x in seen or seen_add(x))]
 
-def uniq2(seq):
-    """ filter out duplicated elements and return the unique ones, slower """
-    res = []
-    for test in seq:
-        if not any(test == c for c in res):
-            res.append(test)
-    return res
-
 def get_partial_derivs(function, variables, n, flatlist=True):
     """ Evaluate the first n partial derivatives of a function with respect to
         specified variables and return a list of unique derivatives """
-    
-    ders = [uniq1(sp.derive_by_array(function, variables))]
+
+    ders = [uniq(sp.derive_by_array(function, variables))]
     for j in range(1, n):
-        ders.append(uniq1(sp.derive_by_array(ders[j-1], variables)))
+        ders.append(uniq(sp.derive_by_array(ders[j-1], variables)))
     return flatten(ders) if flatlist else ders
 
 def get_all_partial_derivs(function, variables, n, flatlist=True):
diff --git a/classical/harmonic.py b/rotagaporp/harmonic.py
similarity index 100%
rename from classical/harmonic.py
rename to rotagaporp/harmonic.py
diff --git a/classical/liouville.py b/rotagaporp/liouville.py
similarity index 100%
rename from classical/liouville.py
rename to rotagaporp/liouville.py
diff --git a/classical/mpliouville.py b/rotagaporp/mpliouville.py
similarity index 100%
rename from classical/mpliouville.py
rename to rotagaporp/mpliouville.py
diff --git a/classical/newtonian.py b/rotagaporp/newtonian.py
similarity index 100%
rename from classical/newtonian.py
rename to rotagaporp/newtonian.py
diff --git a/classical/pbc.py b/rotagaporp/pbc.py
similarity index 100%
rename from classical/pbc.py
rename to rotagaporp/pbc.py
diff --git a/classical/propagators.py b/rotagaporp/propagators.py
similarity index 100%
rename from classical/propagators.py
rename to rotagaporp/propagators.py
diff --git a/classical/symplectic.py b/rotagaporp/symplectic.py
similarity index 100%
rename from classical/symplectic.py
rename to rotagaporp/symplectic.py
diff --git a/classical/systems.py b/rotagaporp/systems.py
similarity index 100%
rename from classical/systems.py
rename to rotagaporp/systems.py
diff --git a/tests/.gitignore b/tests/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..bee8a64b79a99590d5303307144172cfe824fbf7
--- /dev/null
+++ b/tests/.gitignore
@@ -0,0 +1 @@
+__pycache__
diff --git a/classical/mptests.py b/tests/mp_tests.py
similarity index 97%
rename from classical/mptests.py
rename to tests/mp_tests.py
index b57ecc6d9938eb3f09e75d00cce50a49a5db7425..242b28f0b0b4ac379d3a955108ff30e590d37f8c 100644
--- a/classical/mptests.py
+++ b/tests/mp_tests.py
@@ -2,13 +2,14 @@
 import unittest
 import numpy as np
 import sympy as sp
-from units import LJUnitsElement
-from systems import ManyAtomsSystem, MPPS1D, MPPS3D
 from ase import Atoms
-from chebyshev import Chebyshev
-from newtonian import Newtonian
-from symplectic import Symplectic
-from mpliouville import MPLPQP, MPLPAQ, MPLPVQ, MPLPDR, MPLPDV, MPLPNR
+from rotagaporp.systems import ManyAtomsSystem, MPPS1D, MPPS3D
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.newtonian import Newtonian
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.mpliouville import MPLPQP, MPLPAQ, MPLPVQ
+from rotagaporp.mpliouville import MPLPDR, MPLPDV, MPLPNR
+from utilities.units import LJUnitsElement
 
 class MPLPTest(unittest.TestCase):
     """ Comparison tests with different MPLP classes without dynamics """
@@ -46,7 +47,7 @@ class MPLPTest(unittest.TestCase):
 
     def test_qvec_pvec(self):
         """ cross-check all classes with all numerics and variants """
-        from newtonian import leja_points_ga
+        from rotagaporp.newtonian import leja_points_ga
 
         qval = np.random.rand(2, 3)
         pval = np.random.rand(2, 3)
@@ -514,8 +515,8 @@ class TwoParticleXDSystemTest(unittest.TestCase):
 
     def lj_2p_xd_systems_test(self):
         """ collision of two LJ particles with collinear initial momenta """
-        from liouville import LiouvillePowersQP
-        from systems import LennardJones
+        from rotagaporp.liouville import LiouvillePowersQP
+        from rotagaporp.systems import LennardJones
         cases = [
             {'propagator': Chebyshev,
              'params_2p': {'tstart': 0., 'tend': 2., 'tstep': 0.01,
@@ -564,8 +565,8 @@ class TwoParticleXDSystemTest(unittest.TestCase):
 
     def morse_2p_xd_systems_test(self):
         """ Two particles in a Morse potential """
-        from liouville import LiouvillePowersQP
-        from systems import Morse
+        from rotagaporp.liouville import LiouvillePowersQP
+        from rotagaporp.systems import Morse
         cases = [
             {'propagator': Chebyshev,
              'params_2p': {'tstart': 0., 'tend': 10., 'tstep': 0.1,
diff --git a/classical/tests.py b/tests/sp_tests.py
similarity index 96%
rename from classical/tests.py
rename to tests/sp_tests.py
index 8cd32a6adf330354ef2974e1987ba4b5a314fa3f..199fbff991b05614d7efb36979fb9e37048ab5b9 100644
--- a/classical/tests.py
+++ b/tests/sp_tests.py
@@ -2,23 +2,27 @@
 import unittest
 import numpy as np
 import sympy as sp
-from liouville import q, p
-from liouville import LiouvillePowers, LiouvillePowersQP, LiouvillePowersNR
-from symplectic import Symplectic
-from chebyshev import Chebyshev
-from newtonian import Newtonian
-from newtonian import leja_points_ga, leja_points_jb
-from newtonian import chebyshev_points
-from newtonian import newton_scalar
-from systems import Harmonic, Anharmonic, Ballistic, Morse, LennardJones
-from systems import ManyParticlePhysicalSystem
-from forcefield import PairFF, FFHarmonic, FFLennardJones, FFMorse, FFCoulomb
-from mpliouville import MPLPQP
-from pbc import get_distance_matrix_lc as dists1
-from pbc import get_distance_matrix_bc as dists2
-from pbc import wrap_positions as wrap
-from combinatorics import cnk_chebyshev_T1, cnk_chebyshev_T2, cnk_chebyshev_T3
-from combinatorics import cnk_chebyshev_U1, cnk_chebyshev_U2, cnk_chebyshev_U3
+from rotagaporp.liouville import q, p
+from rotagaporp.liouville import LiouvillePowers, LiouvillePowersQP, LiouvillePowersNR
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.newtonian import Newtonian
+from rotagaporp.newtonian import leja_points_ga, leja_points_jb
+from rotagaporp.newtonian import chebyshev_points
+from rotagaporp.newtonian import newton_scalar
+from rotagaporp.systems import Harmonic, Anharmonic, Ballistic
+from rotagaporp.systems import Morse, LennardJones
+from rotagaporp.systems import ManyParticlePhysicalSystem
+from rotagaporp.forcefield import PairFF, FFHarmonic
+from rotagaporp.forcefield import FFLennardJones, FFMorse, FFCoulomb
+from rotagaporp.mpliouville import MPLPQP
+from rotagaporp.pbc import get_distance_matrix_lc as dists1
+from rotagaporp.pbc import get_distance_matrix_bc as dists2
+from rotagaporp.pbc import wrap_positions as wrap
+from rotagaporp.combinatorics import cnk_chebyshev_T1, cnk_chebyshev_T2
+from rotagaporp.combinatorics import cnk_chebyshev_T3
+from rotagaporp.combinatorics import cnk_chebyshev_U1, cnk_chebyshev_U2
+from rotagaporp.combinatorics import cnk_chebyshev_U3
 
 class LiouvillePowersTest(unittest.TestCase):
     """ Test of the one-particle Liouville powers class """
@@ -423,7 +427,7 @@ class DividedDifferenceTest(unittest.TestCase):
 
     def divided_difference_test(self):
         """ compare two algorithms for divided differences """
-        from newtonian import divdiff_1, divdiff_2
+        from rotagaporp.newtonian import divdiff_1, divdiff_2
         nodes = np.arange(-1, 1, 0.1)
         funcvals = lambda x: 2*np.exp(3*x)
         diffs1 = divdiff_1(nodes, funcvals(nodes))
@@ -552,7 +556,7 @@ class ChebyshevAssessmentTest(unittest.TestCase):
 
     def chebyshev_scalar_test(self):
         """ chebyshev expansion of the scalar function exp(z*t) """
-        from chebyshev import chebyshev_scalar
+        from rotagaporp.chebyshev import chebyshev_scalar
         from itertools import product
 
         sample_length = 10
@@ -865,7 +869,7 @@ class CombinatoricsTest(unittest.TestCase):
     def test_newton_cnk_int(self):
         """ test Newton expansion coefficients for integer interpolation
             nodes """
-        from combinatorics import cnk_newton_1, cnk_newton_2
+        from rotagaporp.combinatorics import cnk_newton_1, cnk_newton_2
         lambdas = [1, 2, 3, 4]
         cnk1 = cnk_newton_1(lambdas, 4, dtype=int).tolist()
         cnk2 = cnk_newton_2(lambdas, 4, dtype=int).tolist()
@@ -875,7 +879,7 @@ class CombinatoricsTest(unittest.TestCase):
     def test_newton_cnk_float(self):
         """ test Newton expansion coefficients for real interpolation
             nodes """
-        from combinatorics import cnk_newton_1, cnk_newton_2
+        from rotagaporp.combinatorics import cnk_newton_1, cnk_newton_2
         lambdas = np.random.rand(8)
         cnk1 = cnk_newton_1(lambdas, 6, dtype=float)
         cnk2 = cnk_newton_2(lambdas, 6, dtype=float)
@@ -884,7 +888,7 @@ class CombinatoricsTest(unittest.TestCase):
     def test_newton_cnk_complex(self):
         """ test Newton expansion coefficients for complex interpolation
             nodes """
-        from combinatorics import cnk_newton_1, cnk_newton_2
+        from rotagaporp.combinatorics import cnk_newton_1, cnk_newton_2
         reals = np.random.rand(8)
         imags = np.random.rand(8)
         lambdas = reals + 1j * imags
@@ -931,7 +935,7 @@ class GenericFunctionsTests(unittest.TestCase):
 
     def lambdify_long_test(self):
         """ test the lambdify wrapper used for more than 255 arguments """
-        from generic import lambdify_long
+        from rotagaporp.generic import lambdify_long
         import sympy as sp
 
         dim = 64
diff --git a/classical/test_pbc.py b/tests/test_pbc.py
similarity index 88%
rename from classical/test_pbc.py
rename to tests/test_pbc.py
index 0dba6a30fdacdd95e08d6bf38273d24345b182b5..3df3270ee006636de1880df3495995d355455fe6 100644
--- a/classical/test_pbc.py
+++ b/tests/test_pbc.py
@@ -1,16 +1,18 @@
+""" """
 import unittest
 import numpy as np
-from pbc import get_distance_matrix_lc as dists1
-from pbc import get_distance_matrix_bc as dists2
-from pbc import wrap_positions as wrap1
 from ase.geometry import wrap_positions as wrap2
 from ase.geometry import get_distances
+from rotagaporp.pbc import get_distance_matrix_lc as dists1
+from rotagaporp.pbc import get_distance_matrix_bc as dists2
+from rotagaporp.pbc import wrap_positions as wrap1
 
 class PBCTest(unittest.TestCase):
-
+    """ """
     pbc = (True, True, True)
 
     def test_position_small(self):
+        """ """
         cell = np.array([3.0, 4.0, 5.0])
         posvec = np.array([[3.5, 2.0, 2.5], [0.0, -1.0, 4.0]])
         res1 = wrap1(posvec, cell)
@@ -18,6 +20,7 @@ class PBCTest(unittest.TestCase):
         self.assertTrue(np.allclose(res1, res2))
 
     def test_distance_small(self):
+        """ """
         cell = np.array([3.0, 4.0, 5.0])
         posvec = np.array([[1.5, 2.0, 2.5], [2.7, 1.5, 4.3], [1.2, 0.3, 4.2]])
         res1 = np.linalg.norm(dists1(posvec, cell), axis=2)
@@ -30,6 +33,7 @@ class PBCTest(unittest.TestCase):
         self.assertTrue(np.allclose(res1, res_ase))
 
     def test_distance_big(self):
+        """ """
         cell = np.random.rand(3)
         posvec = np.random.rand(100, 3)
         ful1 = dists1(posvec, cell)
diff --git a/utilities/.gitignore b/utilities/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..bee8a64b79a99590d5303307144172cfe824fbf7
--- /dev/null
+++ b/utilities/.gitignore
@@ -0,0 +1 @@
+__pycache__
diff --git a/utilities/__init__.py b/utilities/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..18bc67f60e0baeaa47e724aa01782cad252e00ec
--- /dev/null
+++ b/utilities/__init__.py
@@ -0,0 +1,3 @@
+import os, sys
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
+
diff --git a/classical/assess.py b/utilities/assess.py
similarity index 92%
rename from classical/assess.py
rename to utilities/assess.py
index 42d0c836004a8547ff52dd73b15c84fe8f1200b6..7f5e714b4501921bd64209759f4adb24180596f9 100644
--- a/classical/assess.py
+++ b/utilities/assess.py
@@ -1,3 +1,4 @@
+""" Functions to assess the performance and accuracy of propagators """
 import time
 import itertools
 import numpy as np
@@ -7,6 +8,7 @@ from uuid import uuid4
 from datetime import datetime
 
 def assess_single(propagator, syst, params, l2error=False):
+    """ Start and evaluate one parameter multiplet """
     import traceback
     res_columns = {}
     try:
@@ -38,6 +40,7 @@ def assess_single(propagator, syst, params, l2error=False):
 def assess_multi(propagator, system, params, var_columns, var_columns_names,
         res_columns_names, con_columns, con_columns_names, indx='time stamp',
         l2error=False):
+    """ Generate multiplets from parameters ranges and start evaluations """
     data_values = []
     data_index = []
     for multiplet in itertools.product(*var_columns):
diff --git a/classical/pandasplot.py b/utilities/pandasplot.py
similarity index 100%
rename from classical/pandasplot.py
rename to utilities/pandasplot.py
diff --git a/utilities/scripts/ase-argon.py b/utilities/scripts/ase-argon.py
new file mode 100644
index 0000000000000000000000000000000000000000..a916531948d806e7e0b6af3c3a574adea5ccca51
--- /dev/null
+++ b/utilities/scripts/ase-argon.py
@@ -0,0 +1,38 @@
+""" Molecular dynamics of argon at constant energy """
+
+from ase import units
+from ase.build import bulk
+from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
+from ase.calculators.lj import LennardJones
+from ase.md.verlet import VelocityVerlet
+from ase.md import MDLogger
+from ase.io.trajectory import Trajectory
+
+T0_lj = 3.0 # initial temperature in Ar LJ units
+rho_lj = 0.88 # density in Ar LJ units
+sigma_a = 3.405*units.Angstrom # in Angstrom for Ar
+epsilon_a = 1.65e-21*units.J # in Joule for Ar
+T = T0_lj * epsilon_a / units.kB # initial temperature in Kelvin
+rho = rho_lj/sigma_a**3 # density in atoms per cubic Angstrom
+n_cell = 3 # number of unit cells
+n_atoms = 4*n_cell**3 # 4 atoms in a unit cell
+l_box = (n_atoms/rho)**(1/3) # edge length of simulation box
+l_cell = l_box / n_cell # latice parameter
+
+atoms = bulk('Ar', 'fcc', a=l_cell, cubic=True).repeat((n_cell, n_cell, n_cell))
+rho_real = sum(atoms.get_masses())/atoms.get_volume()*units.m**3/units.kg
+print('density, kg/m^3', rho_real)
+
+lj_params = {'sigma': sigma_a, 'epsilon': epsilon_a}
+MaxwellBoltzmannDistribution(atoms, T*units.kB)
+atoms.set_calculator(LennardJones(**lj_params))
+
+traj = Trajectory('argon.traj', 'w', atoms)
+dyn = VelocityVerlet(atoms, 1.0*units.fs)
+logger = MDLogger(dyn, atoms, 'argon.log', header=True, stress=True, mode='w')
+dyn.attach(logger, interval=100)
+dyn.attach(traj.write, interval=100)
+
+logger() # log initial state
+traj.write(atoms) # write initial structure
+dyn.run(10000)
diff --git a/utilities/scripts/ase-copper.py b/utilities/scripts/ase-copper.py
new file mode 100644
index 0000000000000000000000000000000000000000..ccec7a608ef3e5cd0568e770af395e8867f056da
--- /dev/null
+++ b/utilities/scripts/ase-copper.py
@@ -0,0 +1,41 @@
+"""Demonstrates molecular dynamics with constant energy."""
+
+from ase.lattice.cubic import FaceCenteredCubic
+from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
+from ase.md.verlet import VelocityVerlet
+from ase import units
+from ase.calculators.emt import EMT
+from ase.visualize import view
+
+size = 3
+
+# Set up a crystal
+atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
+                          symbol='Cu',
+                          size=(size, size, size),
+                          pbc=True)
+
+view(atoms)
+
+# Describe the interatomic interactions with the Effective Medium Theory
+atoms.set_calculator(EMT())
+
+# Set the momenta corresponding to T=300K
+MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
+
+# We want to run MD with constant energy using the VelocityVerlet algorithm.
+dyn = VelocityVerlet(atoms, 5 * units.fs)  # 5 fs time step.
+
+
+def printenergy(a):
+    """Function to print the potential, kinetic and total energy"""
+    epot = a.get_potential_energy() / len(a)
+    ekin = a.get_kinetic_energy() / len(a)
+    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
+          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))
+
+# Now run the dynamics
+printenergy(atoms)
+for i in range(20):
+    dyn.run(10)
+    printenergy(atoms)
diff --git a/classical/assess_chebyshev.py b/utilities/scripts/assess_chebyshev.py
similarity index 100%
rename from classical/assess_chebyshev.py
rename to utilities/scripts/assess_chebyshev.py
diff --git a/utilities/scripts/assess_test.py b/utilities/scripts/assess_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..c221fd214a62a99fca529c9accfc81f587e881ab
--- /dev/null
+++ b/utilities/scripts/assess_test.py
@@ -0,0 +1,67 @@
+""" Effort-Accuracy of Chebyshev propagator for an N-particle 3D system """
+
+import numpy as np
+from mpliouville import MPLPAQ
+from chebyshev import Chebyshev
+from pandas.io.json import json_normalize
+from systems import MPPS3D
+import time
+
+# system parameters
+numeric = {'tool': 'lambdify', 'modules': ['mpmath']}
+# numeric = {'tool': 'subs'}
+masses = np.array([2., 2.])
+qval = np.array([[1., 0., 0.], [4., 0., 0.]])
+pval = np.array([[0., 0., 0.], [0., 0., 0.]])
+symbols = ['0x:z', '1x:z']
+ffparams = [{'class': 'FFMorse', 'pair': ((0, 1, 2), (3, 4, 5)),
+             'params': {'distance': 1.}}]
+syspar = {'fpprec': 'multi', 'prec': 30, 'numeric': numeric}
+
+# propagator parameters
+variant = None
+mplp = MPLPAQ
+propagator = 'chebyshev'
+params = {'tstart': 0., 'tend': 10., 'delta': 1., 'fpprec': 'multi', 'prec': 30}
+tstep = 0.01
+nterms = 2
+
+statistics = []
+stats = {}
+stats['propagator'] = propagator
+stats['class'] = mplp.__name__
+stats['variant'] = variant
+stats['nterms'] = nterms
+stats['tstep'] = tstep
+stats['numeric tool'] = numeric['tool']
+stats['nparticles'] = 2
+stats['delta'] = params['delta']
+cheb_params = params.copy()
+cheb_params['nterms'] = nterms
+cheb_params['tstep'] = tstep
+cheb_params['liouville'] = mplp
+cheb_params['syst'] = MPPS3D(ffparams, symbols, masses, qval, pval, **syspar)
+cheb_params['liou_params'] = {'numeric': numeric}
+ts = time.time()
+cheb = Chebyshev(**cheb_params)
+te = time.time()
+stats['preparation time'] = te - ts
+stats['alphan'] = cheb.efficiency
+ts = time.time()
+cheb.propagate()
+te = time.time()
+cheb.analyse()
+runtime = te - ts
+stats['propagation time'] = runtime
+print('cheb.er', cheb.er)
+absvals = np.fabs(cheb.er)
+print('absvals', absvals)
+stats['energy drift'] = np.max(absvals)
+print('edrift', stats['energy drift'])
+stats['time step efficiency'] = tstep/runtime
+rte = (params['tend']-params['tstart'])/runtime
+stats['run time efficiency'] = rte
+statistics.append(stats)
+
+df = json_normalize(statistics)
+df.to_json('time_cheb_subs_testing_only.json', double_precision=15)
diff --git a/classical/assess_verlet.py b/utilities/scripts/assess_verlet.py
similarity index 100%
rename from classical/assess_verlet.py
rename to utilities/scripts/assess_verlet.py
diff --git a/classical/l2_error.py b/utilities/scripts/l2_error.py
similarity index 100%
rename from classical/l2_error.py
rename to utilities/scripts/l2_error.py
diff --git a/utilities/scripts/pandas_merge.py b/utilities/scripts/pandas_merge.py
new file mode 100644
index 0000000000000000000000000000000000000000..8217e7b1a8538ff87c420b975fd9d02a4cf194a3
--- /dev/null
+++ b/utilities/scripts/pandas_merge.py
@@ -0,0 +1,7 @@
+import sys
+import pandas as pd
+
+dataframes = [pd.read_json(f, precise_float=True) for f in sys.argv[1:]]
+df = pd.concat(dataframes, sort=False, ignore_index=True)
+# df = pd.concat(dataframes, sort=False)
+df.to_json('results_merged.json', double_precision=15)
diff --git a/utilities/scripts/plot_contour.py b/utilities/scripts/plot_contour.py
new file mode 100644
index 0000000000000000000000000000000000000000..dd03577544aa0da67b24b30dbd05f33a467c6ba1
--- /dev/null
+++ b/utilities/scripts/plot_contour.py
@@ -0,0 +1,55 @@
+import matplotlib
+matplotlib.use('Qt5Cairo')
+import matplotlib.pyplot as plt
+from scipy.interpolate import griddata
+import numpy as np
+import pandas as pd
+
+files = [
+    'chebyshev_lj_delta/tstep=0.05/results_merged.json',
+#    'chebyshev_anharmonic_high_delta_highprec/prec=60/results_merged.json',
+#    'chebyshev_morse_high_delta_highprec/prec=60/results_merged.json',
+]
+
+dataframes = [pd.read_json(f, precise_float=True) for f in files]
+df = pd.concat(dataframes, sort=False)
+df = df[np.isfinite(df['energy drift'])]
+df['total efficiency'] = np.log10(df['run time efficiency'])-np.log10(df['energy drift'])
+
+df3 = df[(df['tstep'] == 0.05) &
+         (df['propagator'] == 'chebyshev') &
+         (df['energy drift'] < 1.e0)]
+'''
+df3 = df[(df['total efficiency'] > 2) &
+         (df['delta'] == 12.0) & 
+         (df['system'] == 'Morse') &
+         (df['propagator'] == 'newtonian')]
+'''
+# x = df3['delta']
+x = df3['alphan']
+y = df3['nterms']
+# z = df3['total efficiency']
+# z1 = df3['run time efficiency']
+z2 = np.log10(df3['energy drift'])
+# xi = np.linspace(0, 200, 20)
+xi = np.linspace(0, 0.8, 20)
+yi = np.linspace(10, 160, 20)
+# zi = griddata((x, y), z, (xi[None, :], yi[:, None]), method='cubic')
+# zi1 = griddata((x, y), z1, (xi[None, :], yi[:, None]), method='cubic')
+zi2 = griddata((x, y), z2, (xi[None, :], yi[:, None]), method='cubic')
+# zi2 = griddata((x, y), z2, (xi[None, :], yi[:, None]), method='linear')
+
+fig = plt.figure()
+pltobj = fig.add_subplot(111)
+# contour = pltobj.contour(xi, yi, zi)
+# contour1 = plt.contour(xi, yi, zi1)
+# contour2 = pltobj.contour(xi, yi, zi2, colors='g')
+contour2 = pltobj.contour(xi, yi, zi2)
+# plt.clabel(contour, contour.levels)
+# plt.clabel(contour1, contour1.levels)
+plt.clabel(contour2, contour2.levels, fmt='%1d')
+# pltobj.set_xlabel(r'$\Delta L$')
+pltobj.set_xlabel(r'$\alpha/N$')
+pltobj.set_ylabel(r'$N$')
+# plt.legend()
+plt.show()
diff --git a/utilities/scripts/plot_xy.py b/utilities/scripts/plot_xy.py
new file mode 100644
index 0000000000000000000000000000000000000000..7e355011a160ea18f81ebcd11acf75287ddb4e8c
--- /dev/null
+++ b/utilities/scripts/plot_xy.py
@@ -0,0 +1,145 @@
+from pandasplot import PandasPlot
+import pandas as pd
+
+var_labels = {
+    'totals.operations': 'number of operations',
+    'totals.storage': 'number of objects',
+    'preparation time': 'preparation time',
+    'propagation time': 'propagation time',
+    'chebyshev propagation time': 'propagation time',
+    'class': 'MPLP class',
+    'nparticles': 'N particles',
+    'variant': 'variant',
+    'delta': r'$\Delta L$',
+    'nterms': r'Number of terms $N$',
+    'tstep': r'$\Delta t$',
+    'alphan': r'$\alpha/N$',
+    'time step efficiency': r'$\Delta t/T_{runtime}$',
+    'run time efficiency': r'$(t_{end}-t_{start})/T_{run time}$',
+    'energy drift': r'$max(|\Delta E/ E|)$',
+    'L^2 error': r'$max(|q-q_{ref}|)$',
+    'propagator': 'propagator',
+    'numeric tool': 'numeric tool',
+    'max edrift': r'$\Delta\tilde{E}_{t_{end}}$'
+}
+
+files = [
+#    'chebyshev_morse_high_highprec/results.json',
+#    'chebyshev_morse_high_highprec/l2_error/results.json',
+#    'chebyshev_morse_high_more/results.json',
+#    'newtonian_morse_high_highprec/results.json',
+#    'newtonian_morse_high_more/results.json',
+#    'symplectic_morse_high/results.json',
+#    'symplectic_morse_high_highprec/results.json',
+#    'symplectic_morse_high_highprec/l2_error/results.json',
+#    'chebyshev_anharmonic_high_delta_highprec/prec=30/results_merged.json',
+#    'chebyshev_anharmonic_high_delta_highprec/prec=60/results_merged.json',
+#    'chebyshev_anharmonic_high_delta_highprec/prec=60_reproduce/results_merged.json',
+#    'chebyshev_morse_high_delta_highprec/prec=30/results_merged.json',
+    'chebyshev_morse_high_delta_highprec/prec=60/tstep/results_merged.json',
+#    'chebyshev_morse_high_delta_highprec/delta=100/results_merged.json',
+#    'newtonian_morse_high_delta_highprec/new_sets/results_merged.json',
+#    'chebyshev_lj_delta/tstep=0.01/results_merged.json',
+#    'time_cheb_lj_2p_3d_theano_tstep/results_merged.json',
+#    'time_verlet_lj_2p_3d_theano_tstep/results.json',
+#    'time_cheb_lj_2p_3d_lambdify_tstep_highprec/results_merged.json',
+#    'time_verlet_lj_2p_3d_lambdify_tstep_highprec/results.json',
+#    'time_cheb_lj_np_3d_corr/results_merged.json',
+#    'time_cheb_lj_np_3d_lambdify/results_merged.json',
+#    'time_newt_lj_np_3d_lambdify/results_merged.json',
+#    'time_iln_lj_np_3d_subs/results_merged.json',
+#    'time_cheb_mo_2p_3d_lambdify_tstep_highprec/results_merged.json',
+#    'time_verlet_mo_2p_3d_lambdify_tstep_highprec/results.json'
+#    'chebyshev_morse_high_longtime_highprec/alldata.json'
+]
+
+dataframes = [pd.read_json(f, precise_float=True) for f in files]
+df = pd.concat(dataframes, sort=False)
+
+showx = 'nterms'
+showy = 'energy drift'
+
+# add first set
+# select = {'class': 'MPLPDV', 'nparticles': 2}
+# select = {'propagator': 'chebyshev', 'nterms': 30.}
+select = {'propagator': 'chebyshev', 'tstep': 0.5}
+# select = {'propagator': 'chebyshev', 'nterms': 10}
+# select = {'propagator': 'chebyshev', 'class': 'MPLPAQ'}
+# select = {'propagator': 'chebyshev', 'variant': 'dr', 'nparticles': 5}
+# select = {'propagator': 'newtonian', 'tstep': 0.5}
+# select = {'propagator': 'newtonian', 'delta': 12.}
+# select = {'propagator': 'newtonian'}
+# select = {'propagator': 'symplectic'}
+# select = {'numeric tool': 'theano', 'tstep': 0.00001}
+# select = {'variant': 'dv', 'nparticles': 5}
+# select = {'variant': 'dv', 'nterms': 5}
+# select = {'nparticles': 2}
+# legendvars = ['propagator']
+# legendvars = ['alphan']
+# legendvars = ['tstep']
+legendvars = ['delta']
+# legendvars = ['nterms']
+# legendvars = ['nparticles']
+# legendvars = ['variant']
+# legendvars = ['delta']
+# legendvars = ['numeric tool']
+# legendvars = ['class']
+# legendvars = []
+
+pplt1 = PandasPlot(df, showx=showx, showy=showy, labels=var_labels)
+vary = {v: sorted(list(set(pplt1.df[v].dropna().tolist()))) for v in legendvars}
+# pplt1.add_datasets(select=select, vary=vary, miny=1e-32, maxy=1, maxx=1)
+# pplt1.add_datasets(select=select, vary=vary, minx=1e-3, maxx=10, miny=1e-10, maxy=1e4)
+# pplt1.add_datasets(select=select, vary=vary, minx=1e-5, maxx=1e-2, miny=1e-10, maxy=1e4)
+# pplt1.add_datasets(select=select, vary=vary, minx=2, miny=1e-16, maxy=100.)
+# pplt1.add_datasets(select=select, vary=vary, minx=2)
+pplt1.add_datasets(select=select, vary=vary)
+
+# add second set
+# select = {'propagator': 'chebyshev', 'class': 'MPLPNR', 'variant': 'dv'}
+# select = {'numeric tool': 'lambdify', 'tstep': 0.00001}
+# select = {'propagator': 'chebyshev', 'tstep': 1.}
+# select = {'propagator': 'symplectic'}
+# select = {'propagator': 'verlet'}
+# select = {'variant': 'r', 'nparticles': 5}
+# select = {'variant': 'r', 'nterms': 5}
+
+# legendvars = ['nterms']
+# legendvars = []
+# vary = {v: sorted(list(set(pplt1.df[v].dropna().tolist()))) for v in legendvars}
+# pplt1.add_datasets(select=select, vary=vary, marker='s', minx=1e-3, maxx=10)
+# pplt1.add_datasets(select=select, vary=vary, marker='s', minx=1e-2, maxx=10)
+# pplt1.add_datasets(select=select, vary=vary, marker='s', minx=1e-5, maxx=1e-2)
+# pplt1.add_datasets(select=select, vary=vary, marker='s', minx=2)
+# pplt1.add_datasets(select=select, vary=vary, marker='s')
+
+# select = {'class': 'MPLPAQ', 'nterms': 5}
+# select = {'class': 'MPLPAQ', 'nparticles': 5}
+# pplt1.add_datasets(select=select, vary=vary, minx=2)
+# select = {'variant': 'dr', 'nterms': 5}
+# select = {'variant': 'dr', 'nparticles': 5}
+# select = {'propagator': 'chebyshev', 'variant': 'dv', 'nparticles': 5}
+# pplt1.add_datasets(select=select, vary=vary)
+# select = {'propagator': 'chebyshev', 'variant': 'r', 'nparticles': 5}
+# pplt1.add_datasets(select=select, vary=vary)
+
+pplt2 = pplt1
+
+# add further files
+# df = pd.read_json('symplectic_morse_high_highprec/results.json',
+#                    precise_float=True)
+# pplt2 = PandasPlot(df, showx=showx, showy=showy, labels=var_labels,
+#                    plotobj=pplt1.plotobj, colors=pplt1.colors)
+# select = {'propagator': 'chebyshev', 'nterms': 4}
+# legendvars = ['nterms']
+# legendvars = ['variant']
+# legendvars = ['class']
+# vary = {v: sorted(list(set(pplt2.df[v].dropna().tolist()))) for v in legendvars}
+# pplt2.add_datasets(select=select, vary=vary)
+# select = {'propagator': 'chebyshev', 'nterms': 5, 'class': 'MPLPDV'}
+# pplt2.add_datasets(select=select, vary=vary)
+
+
+pplt2.set_axis_type(xtype='log', ytype='log')
+pplt2.show(legend=True)
+# pplt2.show()
diff --git a/utilities/scripts/time_iterjoins.py b/utilities/scripts/time_iterjoins.py
new file mode 100644
index 0000000000000000000000000000000000000000..395d4c8e54d4e570311b2be2b373855abb53d717
--- /dev/null
+++ b/utilities/scripts/time_iterjoins.py
@@ -0,0 +1,53 @@
+import time
+import numpy as np
+import sympy as sp
+from itertools import chain
+
+# list1 = list(range(1000000))
+# list2 = list(range(2000000))
+
+list1 = sp.Array(sp.symbols(' '.join(['q'+str(s) for s in range(1000000)])))
+list2 = sp.Array(sp.symbols(' '.join(['p'+str(s) for s in range(2000000)])))
+
+list1 = np.array(list1)
+list2 = np.array(list2)
+
+#list1 = np.random.rand(1000000)
+#list2 = np.random.rand(2000000)
+
+timings = {'append': 0.0, 'add': 0.0, 'extend': 0.0, 'unpack': 0.0, 'chain': 0.0, 'concatenate': 0.0}
+
+for repeat in range(1):
+#    start = time.time()
+#    for el in list2:
+#        list1.append(el)
+#    end = time.time()
+#    timings['append'] += (end-start)
+
+#    start = time.time()
+#    list1.extend(list2)
+#    end = time.time()
+#    timings['extend'] += (end-start)
+
+#    start = time.time()
+#    list3 = list1 + list2
+#    end = time.time()
+#    print(list3)
+#    timings['add'] += (end-start)
+
+    start = time.time()
+    list4 = [*list1, *list2]
+    end = time.time()
+    timings['unpack'] += (end-start)
+
+    start = time.time()
+    list5 = list(chain(list1, list2))
+    end = time.time()
+    timings['chain'] += (end-start)
+
+    start = time.time()
+    list6 = np.concatenate((list1, list2))
+    end = time.time()
+    timings['concatenate'] += (end-start)
+
+print(timings)
diff --git a/utilities/scripts/time_uniq.py b/utilities/scripts/time_uniq.py
new file mode 100644
index 0000000000000000000000000000000000000000..778ae68b2ddfab5503495e7b386f84f3fde546d1
--- /dev/null
+++ b/utilities/scripts/time_uniq.py
@@ -0,0 +1,29 @@
+import numpy as np
+from utilities.timing import timeit
+
+def uniq1(seq):
+    """ filter out duplicated elements and return the unique ones, faster """
+    seen = set()
+    seen_add = seen.add
+    return [x for x in seq if not (x in seen or seen_add(x))]
+
+def uniq2(seq):
+    """ filter out duplicated elements and return the unique ones, slower """
+    res = []
+    for test in seq:
+        if not any(test == c for c in res):
+            res.append(test)
+    return res
+
+@timeit
+def time_uniq1(seq):
+    return uniq1(seq)
+
+@timeit
+def time_uniq2(seq):
+    return uniq2(seq)
+
+for irange in [10, 100]:
+    print('integers range', irange)
+    seq = np.random.randint(0, irange, (1000000, ))
+    assert time_uniq1(seq) == time_uniq2(seq)
diff --git a/utilities/tests/test_decorator.py b/utilities/tests/test_decorator.py
new file mode 100644
index 0000000000000000000000000000000000000000..f82e58736e22478c8d10f4c0dc4572c2e5086ded
--- /dev/null
+++ b/utilities/tests/test_decorator.py
@@ -0,0 +1,15 @@
+def def_float(fpprec):
+    def def_float_fpprec(f):
+        if fpprec == 'fixed':
+            def func2(*args, **kwargs):
+                return float(f(*args, **kwargs))
+        else:
+            func2 = f
+        return func2
+    return def_float_fpprec
+
+@def_float('fixed')
+def func(x):
+    return x**2
+    
+print(func(2))
diff --git a/utilities/tests/test_griddata.py b/utilities/tests/test_griddata.py
new file mode 100644
index 0000000000000000000000000000000000000000..697ae76c432d1688ee23e978f824d2732823a7d7
--- /dev/null
+++ b/utilities/tests/test_griddata.py
@@ -0,0 +1,31 @@
+from numpy.random import uniform, seed
+# from matplotlib.mlab import griddata
+from scipy.interpolate import griddata
+import matplotlib.pyplot as plt
+import numpy as np
+# make up data.
+npts = 200
+x = uniform(-2, 2, npts)      # 200 random points between -2 to 2
+y = uniform(-2, 2, npts)
+z = x*np.exp(-x**2 - y**2)
+# define grid.
+xi = np.linspace(-2.1, 2.1, 100)	# 100 x grid points
+yi = np.linspace(-2.1, 2.1, 200)	# 200 y grid points
+
+# mlab.griddata accepts the xi and yi 1d arrays above with different lengths. 
+# scipy's griddata requires a full grid array. This is created with the mgrid function. 
+xi, yi = np.mgrid[-2.1:2.1:100j, -2.1:2.1:200j]
+# grid the data.
+# points = np.vstack((x,y)).T
+zi = griddata((x,y), z, (xi, yi), method='cubic')
+print(zi.tolist())
+# contour the gridded data, plotting dots at the nonuniform data points.
+CS = plt.contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
+CS = plt.contourf(xi, yi, zi, 15)
+plt.colorbar()  # draw colorbar
+# plot data points.
+# plt.scatter(x, y, marker='o', s=5, zorder=10)
+plt.xlim(-2, 2)
+plt.ylim(-2, 2)
+plt.title('griddata test (%d points)' % npts)
+plt.show()
diff --git a/utilities/tests/test_prec.py b/utilities/tests/test_prec.py
new file mode 100644
index 0000000000000000000000000000000000000000..72d72b98f29957d63169da0f03258acd5f1ca11c
--- /dev/null
+++ b/utilities/tests/test_prec.py
@@ -0,0 +1,16 @@
+import numpy as np
+from pandas.io.json import json_normalize
+
+statistics = []
+stats = {}
+en0 = np.array(0.5, dtype='float128')
+ene = 0.5 - 1e-16*np.array(np.random.rand(10000), dtype='float128')
+err = (ene-en0)/en0
+absvals = np.fabs(err)
+maxerr = np.max(absvals)
+print(maxerr)
+stats['maxerr'] = maxerr
+statistics.append(stats)
+
+df = json_normalize(statistics)
+df.to_json('test_prec.json', double_precision=15)
diff --git a/classical/timing.py b/utilities/timing.py
similarity index 82%
rename from classical/timing.py
rename to utilities/timing.py
index f9fe94fff05cf9428543260b35dcd026215c6694..d22c16c2110acd8e716e9edafd4fc3b7cdd92272 100644
--- a/classical/timing.py
+++ b/utilities/timing.py
@@ -1,4 +1,6 @@
+""" A simple timing function that can be used as decorator """
 import time
+
 def timeit(method):
     """ return runtime in seconds of a function """
     def timed(*args, **kw):
diff --git a/classical/trajectory.py b/utilities/trajectory.py
similarity index 73%
rename from classical/trajectory.py
rename to utilities/trajectory.py
index 40bf50f43d59133b28206ddd087724ddff8c3ffc..431a6ffe48cfcf9cf612e4696fe35d56ff0f3a91 100644
--- a/classical/trajectory.py
+++ b/utilities/trajectory.py
@@ -1,9 +1,12 @@
+""" Trajectory utility functions """
 from ase import Atoms
 from ase.io import write
 from units import LJUnitsElement
 
-def write_ase_trajectory(propagator, atoms, filename='trajectory.traj', step=1):
-    """ currently restricted to a single element, no export of time axis """
+def write_ase_trajectory(propagator, atoms, filename='trajectory.traj',
+                         step=1):
+    """ propagator is an object of any rotagaporp.Propagator subclass,
+        currently restricted to a single element, no export of time axis """
     symbols = atoms.get_chemical_symbols()
     assert len(set(symbols)) == 1
     ljunits = LJUnitsElement(symbols[0])
diff --git a/classical/units.py b/utilities/units.py
similarity index 100%
rename from classical/units.py
rename to utilities/units.py