diff --git a/examples/test_cheb_mo_1p_1d_symmetric.py b/examples/test_cheb_mo_1p_1d_symmetric.py
index 17254c926982f765da7e3d945fcd306714a92e44..66538973ccc9a6033be01f2aaa579ec52f8523ed 100644
--- a/examples/test_cheb_mo_1p_1d_symmetric.py
+++ b/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}
diff --git a/examples/test_cheb_mo_2p_3d_symmetric.py b/examples/test_cheb_mo_2p_3d_symmetric.py
new file mode 100644
index 0000000000000000000000000000000000000000..08f6df47ff49e3281768f2a258d6f80fff9dd34e
--- /dev/null
+++ b/examples/test_cheb_mo_2p_3d_symmetric.py
@@ -0,0 +1,105 @@
+""" 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()
diff --git a/rotagaporp/liouville.py b/rotagaporp/liouville.py
index 316a23cc645ca8c94e69fd42192210de14c0cf4c..05cf443ca6f6397049134617afaa468df5c24f48 100644
--- a/rotagaporp/liouville.py
+++ b/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 """
diff --git a/rotagaporp/propagators.py b/rotagaporp/propagators.py
index 47184261a1c9a9a9f1fbda7fc65e7b560ad09388..95820e3d9948b55ce7cb06605c737348fc9346c2 100644
--- a/rotagaporp/propagators.py
+++ b/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)
diff --git a/tests/sy_tests.py b/tests/sy_tests.py
index cac8387831d698204fdcbe65c5b76f0a28cf524e..9c09777c4d9786b7a2c253429c9626b8534ecaef 100644
--- a/tests/sy_tests.py
+++ b/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))