self-starting symmetric propagator for many particles, example and test

examples/test_cheb_mo_1p_1d_symmetric.py

@@ -14,7 +14,7 @@ numeric = {'tool': 'lambdify', 'modules': 'mpmath'}
 
 params = {'tstart': spfloat(0., prec), 'tend': spfloat(100., prec),
           'tstep': spfloat(0.2, prec), 'nterms': 6, 'tqdm': True}
-# params = {'tstart': 0., 'tend': 100., 'tstep': 0.4, 'nterms': 4, 'tqdm': True}
+# params = {'tstart': 0., 'tend': 100., 'tstep': 0.2, 'nterms': 6, 'tqdm': True}
 
 syspar = {'q0': spfloat(3., prec), 'p0': spfloat(0., prec),
           'fpprec': 'multi', 'prec': prec}

examples/test_cheb_mo_2p_3d_symmetric.py

+""" Three dimensional Lennard-Jones gas with two particles """
+import numpy as np
+import pickle as pkl
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import SymmetricChebyshev
+from rotagaporp.chebyshev import Chebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.mpliouville import MPLPQP
+from rotagaporp.systems import MPPS3D, Morse
+
+M = [2.0, 2.0]
+numeric = {'tool': 'lambdify', 'modules': 'numpy'}
+params = {'tstart': 0.0, 'tend': 1000., 'tstep': 0.1, 'tqdm': False}
+syspar = {'symbols': ['0x:z', '1x:z'], 'q0': [[1, 2, 3], [4, 2, 3]],
+          'p0': [[0, 0, 0], [0, 0, 0]], 'masses': M, 'numeric': numeric}
+ffpars = [{'class': 'FFMorse', 'pair': ((0, 1, 2), (3, 4, 5)),
+           'params': {'distance': 1.}}]
+system = MPPS3D(ffpars, **syspar)
+
+# velocity Verlet
+vv_prop = Symplectic(syst=system, **params)
+vv_prop.propagate()
+time_vv, q_vv, p_vv = vv_prop.get_trajectory_3d()
+r_vv = np.abs(q_vv[:, 1]-q_vv[:, 0])
+r_vv_1p = np.linalg.norm(r_vv, axis=1)
+pr_vv = (p_vv[:, 1]*M[0]-p_vv[:, 0]*M[1])/np.sum(M)
+pr_vv_1p = np.einsum('ij,ij->i', pr_vv, r_vv/r_vv_1p[:, np.newaxis])
+vv_prop.analyse()
+print('energy drift of velocity verlet prop', np.max(np.fabs(vv_prop.er)))
+
+system_1d = Morse(q0=3., p0=0.)
+q_an, p_an = system_1d.solution(times=time_vv)
+en0 = system_1d.energy(system_1d.q0, system_1d.p0)
+evecf = np.vectorize(system_1d.energy)
+err_an = np.array((evecf(q_an, p_an)-en0)/en0, dtype=float)
+print('energy drift of analytic solution', np.max(np.fabs(err_an)))
+
+# one-step propagator
+params['liouville'] = MPLPQP
+params['nterms'] = 4
+ch_prop = Chebyshev(syst=system, **params)
+ch_prop.propagate()
+time_ch, q_ch, p_ch = ch_prop.get_trajectory_3d()
+r_ch = np.abs(q_ch[:, 1]-q_ch[:, 0])
+r_ch_1p = np.linalg.norm(r_ch, axis=1)
+pr_ch = (p_ch[:, 1]*M[0]-p_ch[:, 0]*M[1])/np.sum(M)
+pr_ch_1p = np.einsum('ij,ij->i', pr_ch, r_ch/r_ch_1p[:, np.newaxis])
+ch_prop.analyse()
+print('energy drift of one-step prop', np.max(np.fabs(ch_prop.er)))
+
+# symmetrized propagator
+params['numeric'] = numeric
+ts_prop = SymmetricChebyshev(syst=system, **params)
+ts_prop.propagate()
+time_ts, q_ts, p_ts = ts_prop.get_trajectory_3d()
+r_ts = np.abs(q_ts[:, 1]-q_ts[:, 0])
+r_ts_1p = np.linalg.norm(r_ts, axis=1)
+pr_ts = (p_ts[:, 1]*M[0]-p_ts[:, 0]*M[1])/np.sum(M)
+pr_ts_1p = np.einsum('ij,ij->i', pr_ts, r_ts/r_ts_1p[:, np.newaxis])
+ts_prop.analyse()
+print('energy drift of symmetric prop', np.max(np.fabs(ts_prop.er)))
+
+assert len(time_ts) == len(time_ch) == len(time_vv)
+assert np.allclose(time_ts, time_ch)
+assert np.allclose(time_ts, time_vv)
+
+fig = plt.figure()
+plot1 = fig.add_subplot(311)
+plot1.plot(time_ch, np.maximum.accumulate(abs(ch_prop.er)),
+           label='Chebyshev')
+plot1.plot(time_ts, np.maximum.accumulate(abs(ts_prop.er)),
+           label='Symmetric Chebyshev')
+plot1.plot(time_ts, np.maximum.accumulate(abs(vv_prop.er)),
+           label='Velocity Verlet')
+plot1.set_ylabel('maximum energy drift')
+plot1.set_yscale('log')
+plot1.set_xscale('log')
+plot1.legend()
+
+plot2 = fig.add_subplot(312)
+plot2.plot(time_ts, r_ch_1p, label='Chebyshev')
+plot2.plot(time_ts, r_ts_1p, label='Symmetric Chebyshev')
+plot2.plot(time_ts, r_vv_1p, label='Velocity Verlet')
+plot2.plot(time_ts, q_an, label='Analytic solution')
+plot2.set_xlabel(r'$t$')
+plot2.set_ylabel(r'$r(t)$')
+plot2.set_xscale('log')
+plot2.legend()
+
+plot3 = fig.add_subplot(313)
+plot3.plot(time_ts, pr_ch_1p, label='Chebyshev')
+plot3.plot(time_ts, pr_ts_1p, label='Symmetric Chebyshev')
+plot3.plot(time_ts, pr_vv_1p, label='Velocity Verlet')
+plot3.plot(time_ts, p_an, label='Analytic solution')
+plot3.set_xlabel(r'$t$')
+plot3.set_ylabel(r'$p_r(t)$')
+plot3.set_xscale('log')
+plot3.legend()
+
+with open('pltobj.pkl', 'wb') as fh:
+    pkl.dump(fig, fh)
+
+fig.savefig('pltfig.png')
+
+# plt.show()

rotagaporp/liouville.py

@@ -218,6 +218,8 @@ class SPLPQP:
                          function=system.q, **kwargs)
         self.pobj = SPLP(nterms, system.hamiltonian, system.q, system.p,
                          function=system.p, **kwargs)
+        self.ilnqcart = self.qobj.Lf
+        self.ilnpcart = self.pobj.Lf
 
     def get_val_float(self, qval, pval):
         """ evaluate the symbolic expressions for given values of q and p """

rotagaporp/propagators.py

@@ -3,6 +3,7 @@
 import mpmath as mp
 import numpy as np
 import sympy as sp
+from collections.abc import Iterable
 from liouville import SPLPQP
 
 class Propagator:
@@ -206,12 +207,17 @@ class SymmetricPolynomialPropagator(PolynomialPropagator):
         multi   mp.Float    subs        None        sympy       step_sympy
         multi   mp.Float    lambdify    mpmath      mpmath      step_mpmath
         -------------------------------------------------------------------
+        NOTABENE: for numeric tool lambdify modules 'tensorflow' and 'math'
+        are currently not supported.
         """
 
         from generic import validate_numeric
 
         if numeric is not None:
             self.numeric = numeric
+        elif hasattr(self.syst, 'numeric'):
+            self.numeric = self.syst.numeric
+
         tool = self.numeric['tool']
         assert tool in ['subs', 'lambdify', 'theano']
         if tool == 'subs':
@@ -247,72 +253,81 @@ class SymmetricPolynomialPropagator(PolynomialPropagator):
             self.matinv = np.linalg.inv
             self.normf = np.linalg.norm
 
+        if isinstance(self.q, Iterable):
+            assert isinstance(self.p, Iterable)
+            assert isinstance(self.syst.q, Iterable)
+            assert isinstance(self.syst.p, Iterable)
+            assert len(self.q) == len(self.p)
+            assert len(self.q) == len(self.syst.q)
+            assert len(self.syst.q) == len(self.syst.p)
+            self.concat = lambda a, b: (*a, *b)
+            self.split = lambda x: (x[:len(x)//2], x[len(x)//2:])
+        else:
+            self.concat = lambda a, b: (a, b)
+            self.split = lambda x: (x[0], x[1])
+
     def set_gvar_jacobian(self):
         """ initialize the G_var vector function and the Jacobian matrix """
 
         if self.fpprec == 'multi':
             self.coeffs = sp.Matrix(self.coeffs)
             self.reverse_coeffs = sp.Matrix(self.reverse_coeffs)
-
-        gq = self.syst.q + self.reverse_coeffs.dot(self.iln.qobj.Lf)
-        gp = self.syst.p + self.reverse_coeffs.dot(self.iln.pobj.Lf)
-        self.z = sp.Matrix((self.syst.q, self.syst.p))
-        jacobian = sp.Matrix([gq, gp]).jacobian(self.z)
+        ugq =self.reverse_coeffs.dot(self.iln.ilnqcart)
+        ugp = self.reverse_coeffs.dot(self.iln.ilnpcart)
+        if isinstance(self.syst.q, Iterable):
+            gq = sp.Matrix(self.syst.q) + sp.Matrix(ugq)
+            gp = sp.Matrix(self.syst.p) + sp.Matrix(ugp)
+        else:
+            gq = self.syst.q + ugq
+            gp = self.syst.p + ugp
+        self.z = sp.Matrix(self.concat(self.syst.q, self.syst.p))
+        jacobian = sp.Matrix(self.concat(gq, gp)).jacobian(self.z)
         self.matrix_shape = jacobian.shape
+        self.vector_shape = (self.z.shape[0],)
         if self.numeric['tool'] == 'lambdify':
             modules = self.numeric['modules']
             if self.numeric['modules'] == 'mpmath':
-                gvar = sp.Matrix([gq, gp])
+                gvar = sp.Matrix(self.concat(gq, gp))
             else:
-                gvar = sp.Array([gq, gp])
+                gvar = sp.Array(self.concat(gq, gp))
             self.jacobian = sp.lambdify(self.z, jacobian, modules)
             self.gvar = sp.lambdify(self.z, gvar, modules)
         elif self.numeric['tool'] == 'theano':
             from sympy.printing.theanocode import theano_function as tf
-            gvar = sp.Matrix([gq, gp])
+            gvar = sp.Matrix(self.concat(gq, gp))
             self.jacobian = tf(self.z, jacobian, on_unused_input='ignore')
             self.gvar = tf(self.z, gvar, on_unused_input='ignore')
         else:
             self.jacobian = jacobian
-            self.gvar = sp.Matrix([gq, gp])
+            self.gvar = sp.Matrix(self.concat(gq, gp))
 
-    def matrix_eval(self, matrix, zval):
+    def array_eval(self, array, shape, zval):
+        """ evaluate the elements of a matrix that are sympy expressions """
         if self.numeric['tool'] == 'subs':
-            ret = matrix.xreplace(dict(zip(self.z, zval))).evalf()
+            retval = array.xreplace(dict(zip(self.z, zval))).evalf()
         else:
+            retval = array(*zval)
             if self.linalg == 'numpy':
-                ret = np.array(matrix(*zval))
-                ret.shape = self.matrix_shape
-            else:
-                ret = matrix(*zval)
-        return ret
-
-    def vector_eval(self, vector, zval):
-        if self.numeric['tool'] == 'subs':
-            ret = vector.xreplace(dict(zip(self.z, zval))).evalf()
-        else:
-            if self.linalg == 'numpy':
-                ret = np.array(vector(*zval))
-            else:
-                ret = vector(*zval)
-        return ret
+                retval = np.array(retval)
+                retval.shape = shape
+        return retval
 
     def step_nopbc(self, time, qval, pval):
-        """ implicit stepping procedure using Newton iteration """
-        zval = self.matrix([qval, pval])
+        """ implicit time stepping procedure using Newton iteration """
         qvec, pvec = self.iln.get_val_float(qval, pval)
-        zval0 = self.matrix([self.coeffs.dot(qvec), self.coeffs.dot(pvec)])
-        gcon = zval + zval0
-        zval1 = zval0.copy()
+        qt = self.coeffs.dot(qvec)
+        pt = self.coeffs.dot(pvec)
+        zval = self.matrix(self.concat(qt, pt))
+        gcon = self.matrix(self.concat(qval, pval)) + zval
         residue = 1.
         numiter = 0
         while residue > 10**(-self.prec/2):
             if numiter > self.maxiter:
                 raise RuntimeError('maximum number of iterations exceeded')
-            jacval = self.matrix_eval(self.jacobian, zval1)
-            gvar = self.vector_eval(self.gvar, zval1)
+            jacval = self.array_eval(self.jacobian, self.matrix_shape, zval)
+            gvar = self.array_eval(self.gvar, self.vector_shape, zval)
             correction = self.matmul(self.matinv(jacval), gvar-gcon)
-            zval1 = zval1 - correction
+            zval = zval - correction
             residue = self.normf(correction)
             numiter += 1
-        return (zval1[0], zval1[1])
+        return self.split(zval)

tests/sy_tests.py

@@ -4,7 +4,8 @@ import numpy as np
 from rotagaporp.generic import spfloat
 from rotagaporp.chebyshev import SymmetricChebyshev
 from rotagaporp.newtonian import SymmetricNewtonian
-from rotagaporp.systems import Morse
+from rotagaporp.systems import Morse, MPPSMD
+from rotagaporp.mpliouville import MPLPQP
 
 class TestSymmetricPolynomialPropagator(unittest.TestCase):
     """ Test the self-starting symmetric polynomial propagator """
@@ -18,7 +19,6 @@ class TestSymmetricPolynomialPropagator(unittest.TestCase):
         qan, pan = syst.solution(times=tsy)
         qanfl = np.array(qan, dtype=float)
         panfl = np.array(pan, dtype=float)
-        print(qsy, qanfl)
         self.assertTrue(np.allclose(qsy, qanfl, atol=atol))
         self.assertTrue(np.allclose(psy, panfl, atol=atol))
         self.assertTrue(np.allclose(sy.er, 0, atol=atol))
@@ -148,3 +148,38 @@ class TestSymmetricPolynomialPropagator(unittest.TestCase):
                   'fpprec': 'multi', 'prec': prec}
         system = Morse(**syspar)
         self.compare(SymmetricNewtonian, params, system, 1e-2)
+
+
+class TestSymmetricPolynomialPropagatorMP(unittest.TestCase):
+    """ Self-starting symmetric polynomial propagator - many particles """
+
+    def test_cheb_morse_2p_2d_lambdify(self):
+        """ Chebyshev for two-particle 2D Morse oscillator with lambdify """
+
+        M = [2.0, 2.0]
+        numeric = {'tool': 'lambdify', 'modules': 'numpy'}
+        params = {'tstart': 0.0, 'tend': 5., 'tstep': 0.1, 'nterms': 4,
+                  'liouville': MPLPQP, 'numeric': numeric}
+        syspar = {'symbols': ['0x:y', '1x:y'], 'q0': [[1, 2], [4, 2]],
+                  'p0': [[0, 0], [0, 0]], 'masses': M, 'numeric': numeric}
+        ffpars = [{'class': 'FFMorse', 'pair': ((0, 1), (2, 3)),
+                   'params': {'distance': 1.}}]
+
+        system = MPPSMD(ffpars, dim=2, **syspar)
+        ts_prop = SymmetricChebyshev(syst=system, **params)
+        ts_prop.propagate()
+        time_ts, q_ts, p_ts = ts_prop.get_trajectory()
+        q_ts.shape = (len(q_ts), -1, 2)
+        p_ts.shape = (len(p_ts), -1, 2)
+        r_ts = np.abs(q_ts[:, 1]-q_ts[:, 0])
+        r_ts_1p = np.linalg.norm(r_ts, axis=1)
+        pr_ts = (p_ts[:, 1]*M[0]-p_ts[:, 0]*M[1])/np.sum(M)
+        pr_ts_1p = np.einsum('ij,ij->i', pr_ts, r_ts/r_ts_1p[:, np.newaxis])
+        ts_prop.analyse()
+
+        system_1d = Morse(q0=3., p0=0.)
+        q_an, p_an = system_1d.solution(times=time_ts)
+
+        self.assertTrue(np.allclose(r_ts_1p, q_an, atol=1e-2))
+        self.assertTrue(np.allclose(pr_ts_1p, p_an, atol=1e-2))
+        self.assertTrue(np.allclose(ts_prop.er, 0, atol=1e-2))