Skip to content
Snippets Groups Projects
Commit e6b8dad4 authored by Ivan Kondov's avatar Ivan Kondov
Browse files

enabled full high precision using spfloat() function and adaptions

parent 71b35f4b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
classical/chebyshev.py
+ 4
− 3
View file @ e6b8dad4
""" Chebyshev propagator for classical particle dynamics """
import numpy as np
import sympy as sp
from scipy.special import jn, iv
from propagators import PolynomialPropagator
from liouville import LiouvillePowersQP
......@@ -14,7 +15,6 @@ def besf(fpprec, signed, prec=None):
if fpprec == 'fixed':
bessel = iv if signed else jn
elif fpprec == 'multi':
import sympy as sp
fargs = [] if prec is None else [prec]
bess = sp.besseli if signed else sp.besselj
def bessel(orders, args):
......@@ -42,16 +42,17 @@ class Chebyshev(PolynomialPropagator):
self.shift = -(2.*self.lambda_min/self.delta+1.)
bf = besf(fpprec, signed, prec=prec)
sign = -1 if signed else 1
exp = np.exp if fpprec == 'fixed' else lambda x: sp.exp(x).evalf(prec)
if self.substep is not None:
self.cheb_coef = []
for stime in self.subtime:
alpha = self.delta*stime/2.
scaling = np.exp((self.lambda_min+self.delta/2.0)*stime)
scaling = exp((self.lambda_min+self.delta/2.0)*stime)
coef = cheb_coef(bf, self.nterms, alpha)*scaling
self.cheb_coef.append(coef)
else:
scaling = np.exp((self.lambda_min+self.delta/2.0)*self.tstep)
scaling = exp((self.lambda_min+self.delta/2.0)*self.tstep)
self.cheb_coef = cheb_coef(bf, self.nterms, alpha)*scaling
self.iln = liouville(self.nterms, self.syst, scale=self.scale,
......
classical/generic.py
+ 10
− 0
View file @ e6b8dad4
......@@ -2,6 +2,7 @@
import sympy as sp
def _float(fpprec, f):
""" function wrapper for conditional casting to python float """
if fpprec == 'fixed':
def newfunc(*args, **kwargs):
return float(f(*args, **kwargs))
......@@ -9,6 +10,15 @@ def _float(fpprec, f):
newfunc = f
return newfunc
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
def flatten(x):
""" Return a flat list from a deeply nested iterables """
from collections import Iterable
......
classical/newtonian.py
+ 8
− 7
View file @ e6b8dad4
""" Newtonian propagator for classical particle dynamics """
import numpy as np
import sympy as sp
from generic import spfloat
from propagators import PolynomialPropagator
from liouville import LiouvillePowersQP
......@@ -20,7 +21,7 @@ def divdiff_2(x, y, dtype=float, fpprec='fixed', prec=None):
""" compute divided differences for y_i=f(x_i), alternative algorithm """
n = len(x)
assert len(x) == len(y)
zero = {'fixed': dtype(0), 'multi': sp.N(dtype(0), prec)}
zero = {'fixed': dtype(0), 'multi': spfloat(0., prec)}
a = np.full((n, n), zero[fpprec])
a[:, 0] = y
for i in range(1, n):
......@@ -32,8 +33,8 @@ def divdiff_2(x, y, dtype=float, fpprec='fixed', prec=None):
def trial_points(number, dtype=float, fpprec='fixed', prec=None):
""" construct the trial points for the algorithm from
G. Ashkenazi et al. JCP 103, 10005 (1995) """
zero = {'fixed': dtype(0), 'multi': sp.N(dtype(0), prec)}
one = {'fixed': dtype(1), 'multi': sp.N(dtype(1), prec)}
zero = {'fixed': dtype(0), 'multi': spfloat(0., prec)}
one = {'fixed': dtype(1), 'multi': spfloat(1., prec)}
ci = {'fixed': 1j, 'multi': sp.N(1j, prec)}
pi = {'fixed': np.pi, 'multi': sp.N(sp.pi, prec)}
cosf = {'fixed': np.cos, 'multi': lambda x: sp.N(sp.cos(x), prec)}
......@@ -65,7 +66,7 @@ def leja_points_ga(number, dtype=float, fpprec='fixed', prec=None):
G. Ashkenazi et al. JCP 103, 10005 (1995) """
isclose = {'fixed': np.isclose,
'multi': lambda x, y: abs(x-y) <= sp.N(10**(-prec/2), prec)}
'multi': lambda x, y: abs(x-y) < 10**(-prec/2)}
ntrials = {float: 3*number, complex: 2*number}
trials, center = trial_points(ntrials[dtype], dtype, fpprec, prec)
......@@ -86,9 +87,9 @@ def leja_points_jb(nflp, dtype=float, fpprec='fixed', prec=None):
assert dtype == float
ndtype = {'fixed': dtype, 'multi': np.object}
scaling = {'fixed': float(1), 'multi': sp.N(float(1), prec)}
left = {'fixed': dtype(-1), 'multi': sp.N(dtype(-1), prec)}
right = {'fixed': dtype(1), 'multi': sp.N(dtype(1), prec)}
scaling = {'fixed': dtype(1), 'multi': spfloat(1., prec)}
left = {'fixed': dtype(-1), 'multi': spfloat(-1., prec)}
right = {'fixed': dtype(1), 'multi': spfloat(1., prec)}
zs = np.empty(nflp, dtype=ndtype[fpprec])
zt = np.empty(nflp, dtype=ndtype[fpprec])
......
classical/propagators.py
+ 6
− 5
View file @ e6b8dad4
......@@ -18,7 +18,8 @@ class Propagator:
self.en0 = self.syst.energy(self.q0, self.p0)
assert self.tstart < self.tend
nsteps = int(np.floor((self.tend-self.tstart)/self.tstep))+1
self.time = np.linspace(self.tstart, self.tend, nsteps)
# self.time = np.linspace(self.tstart, self.tend, nsteps)
self.time = [tstart+(tend-tstart)/(nsteps-1)*s for s in range(nsteps)]
self.sq = [self.q]
self.sp = [self.p]
self.pbc = pbc
......@@ -55,11 +56,11 @@ class Propagator:
""" perform a time step - stub method """
return (coordinate, momentum)
def get_trajectory(self, step=1):
def get_trajectory(self, step=1, dtype=np.float):
""" return a trajectory - a time series of coordinate and momentum """
time = np.array(self.time, dtype='float')[::step]
coordinate = np.array(self.sq, dtype='float')[::step]
momentum = np.array(self.sp, dtype='float')[::step]
time = np.array(self.time, dtype=dtype)[::step]
coordinate = np.array(self.sq, dtype=dtype)[::step]
momentum = np.array(self.sp, dtype=dtype)[::step]
return time, coordinate, momentum
def get_trajectory_3d(self, step=1):
......
classical/systems.py
+ 55
− 52
View file @ e6b8dad4
......@@ -5,25 +5,19 @@ from forcefield import FFList
from pbc import get_distance_matrix_lc as get_distance_matrix
from pbc import wrap_positions
from generic import lambdify_long as lambdify
from generic import _float
from generic import _float, spfloat
class PhysicalSystem:
analytic = False
""" base class for one-particle systems """
def __init__(self, potential_energy, mass, q=sp.Symbol('q'),
p=sp.Symbol('p'), q0=None, p0=None, fpprec='fixed',
prec=None):
self.q = q
self.p = p
if fpprec == 'fixed':
self.q0 = q0
self.p0 = p0
self.mass = mass
elif fpprec == 'multi':
self.q0 = sp.Float(q0, prec)
self.p0 = sp.Float(p0, prec)
self.mass = sp.Float(mass, prec)
else:
raise ValueError('invalid parameter: fpprec == '+fpprec)
assert fpprec in ['fixed', 'multi']
self.q0 = q0 if fpprec == 'fixed' else spfloat(q0, prec)
self.p0 = p0 if fpprec == 'fixed' else spfloat(p0, prec)
self.mass = mass if fpprec == 'fixed' else spfloat(mass, prec)
self.fpprec = fpprec
self.prec = prec
self.kinetic_energy = self.p**2/(2*self.mass)
......@@ -48,46 +42,55 @@ class PhysicalSystem:
def _get_kinetic_energy(self, p):
return self.kinetic_energy.subs(self.p, p)
def solution(self, q0, p0, time):
pass
class Harmonic(PhysicalSystem):
delta = 1.0
omega = 1.0
lambda_min = -delta/2.0
analytic = True
def __init__(self, mass=1, **kwargs):
""" harmonic oscillator """
omega = 1.
def __init__(self, mass=1., **kwargs):
q = sp.Symbol('q')
k = self.omega**2*mass
potential_energy = k*q**2/2
PhysicalSystem.__init__(self, potential_energy, mass=mass, **kwargs)
super().__init__(potential_energy, mass=mass, **kwargs)
def solution(self, q0, p0, time):
v0 = p0/self.mass
A = np.sqrt((self.omega**2*q0**2+v0**2)/self.omega**2)
phi = np.arccos(q0/A)
q_ref = A*np.cos(self.omega*time+phi)
p_ref = -A*self.omega*np.sin(self.omega*time+phi)*self.mass
return (q_ref, p_ref)
def solution(self, q0=None, p0=None, times=None):
""" analytical solution of the equation of motion """
if self.fpprec == 'fixed':
from numpy import sin, cos
from numpy import sqrt, arccos
times = np.array(times, dtype=np.float)
else:
from sympy import sqrt
from sympy import acos as arccos
sin = lambda x: np.array(sp.Array(x).applyfunc(sp.sin))
cos = lambda x: np.array(sp.Array(x).applyfunc(sp.cos))
times = np.array(times, dtype=np.object)
qini = self.q0 if q0 is None else q0
pini = self.p0 if p0 is None else p0
v0 = pini/self.mass
A = sqrt((self.omega**2*qini**2+v0**2)/self.omega**2)
phi = arccos(qini/A)
q_ref = A*cos(self.omega*times+phi)
p_ref = -A*self.omega*sin(self.omega*times+phi)*self.mass
return q_ref, p_ref
class Anharmonic(PhysicalSystem):
delta = 1.0
""" one particle in a double-well potential """
def __init__(self, **kwargs):
q = sp.Symbol('q')
potential_energy = (1/2)*q**4 - q**2
PhysicalSystem.__init__(self, potential_energy, mass=1, **kwargs)
super().__init__(potential_energy, mass=1., **kwargs)
class Ballistic(PhysicalSystem):
delta = 1.0
analytic = True
def __init__(self, force=1.0, **kwargs):
""" one particle in a potential with constant force """
def __init__(self, force=1., **kwargs):
self.cf = force
q = sp.Symbol('q')
potential_energy = self.cf*q
PhysicalSystem.__init__(self, potential_energy, mass=1, **kwargs)
super().__init__(potential_energy, mass=1., **kwargs)
def solution(self, q0, p0, time):
q_ref = -0.5*self.cf*time*time + p0*time + q0
......@@ -96,24 +99,24 @@ class Ballistic(PhysicalSystem):
class LennardJones(PhysicalSystem):
analytic = False
def __init__(self, sigma=1., epsilon=1., mass=1, **kwargs):
""" one particle in Lennard-Jones potential """
def __init__(self, sigma=1., epsilon=1., mass=1., **kwargs):
q = sp.Symbol('q')
potential_energy = 4*epsilon*((sigma/q)**12-(sigma/q)**6)
PhysicalSystem.__init__(self, potential_energy, mass, **kwargs)
super().__init__(potential_energy, mass, **kwargs)
class Morse(PhysicalSystem):
delta = 1.
lambda_min = -1.0
analytic = True
def __init__(self, qmin=1.0, kappa=1.0, d=1.0, mass=1, **kwargs):
self.qmin = qmin
self.kappa = kappa
self.d = d
""" one particle in Morse potential """
def __init__(self, qmin=1., kappa=1., d=1., mass=1., fpprec='fixed',
prec=None, **kwargs):
self.qmin = qmin if fpprec == 'fixed' else spfloat(qmin, prec)
self.kappa = kappa if fpprec == 'fixed' else spfloat(kappa, prec)
self.d = d if fpprec == 'fixed' else spfloat(d, prec)
q = sp.Symbol('q')
expr = d*(sp.exp(-2.0*kappa*(q-qmin)) - 2.0*sp.exp(-kappa*(q-qmin)))
PhysicalSystem.__init__(self, expr, mass, **kwargs)
expr = d*(sp.exp(-2*self.kappa*(q-self.qmin))
-2*sp.exp(-self.kappa*(q-self.qmin)))
super().__init__(expr, mass, fpprec=fpprec, prec=prec, **kwargs)
def solution(self, q0=None, p0=None, times=None):
""" Analytic solution as published in Am. J. Phys. 46, 733 (1978) and
......@@ -125,7 +128,7 @@ class Morse(PhysicalSystem):
from numpy import sin, cos, sinh, cosh
from numpy import arccos, arccosh
isclose = np.isclose
times = np.array(times)
times = np.array(times, dtype=np.float)
else:
from sympy import exp, sqrt, sign
from sympy import acos as arccos
......@@ -135,8 +138,8 @@ class Morse(PhysicalSystem):
log = lambda x: np.array(sp.Array(x).applyfunc(sp.log))
sinh = lambda x: np.array(sp.Array(x).applyfunc(sp.sinh))
cosh = lambda x: np.array(sp.Array(x).applyfunc(sp.cosh))
isclose = lambda x, y: abs(x-y) < sp.N(10**(-self.prec/2), self.prec)
times = np.array([sp.N(t, self.prec) for t in times])
isclose = lambda x, y: abs(x-y) < 10**(-self.prec/2)
times = np.array(times, dtype=np.object)
qini = self.q0 if q0 is None else q0
pini = self.p0 if p0 is None else p0
......@@ -176,8 +179,8 @@ class Morse(PhysicalSystem):
assert sqrtarg > 0 or abs(sqrtarg) < np.finfo(float).resolution
t0 = times[0] + direction*2./gamma*arccosh(sqrt(abs(sqrtarg)))
termq = log((sqrt(rho)*cosh(gamma*(times-t0))-1.)/(rho-1.))
termp = (gamma*sqrt(rho)*sinh(gamma*(times-t0))/(sqrt(rho)
*cosh(gamma*(times-t0))-1.))
termp = (gamma*sqrt(rho)*sinh(gamma*(times-t0))
/(sqrt(rho)*cosh(gamma*(times-t0))-1.))
q = self.qmin + termq/self.kappa
p = self.mass*termp/self.kappa
return q, p
......@@ -205,7 +208,7 @@ class ManyParticlePhysicalSystem:
self.p0 = np.array(p0, dtype=dty)
self.mass = np.array(mass, dtype=dty)
elif fpprec == 'multi':
fpfunc = lambda x: sp.Float(x, prec)
fpfunc = lambda x: spfloat(x, prec)
npfunc = np.vectorize(fpfunc, otypes=[np.object])
self.q0 = npfunc(q0)
self.p0 = npfunc(p0)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Help | Imprint | Privacy policy | Accessibility | Contact