from functools import wraps
from numpy import floor, ones_like, einsum, sqrt, absolute, maximum, empty
from numpy.core.multiarray import c_einsum
from mango.debug import debug
from mango.constants import c
[docs]class periodic():
"""
Periodic Boundary Conditions.
The positions are only wrapped at the end of a single step
because the values are not used until they are saved to file.
The distances are always wrapped because they are used during a single step.
"""
__slots__ = ['boxsize', 'pbc', 'posones', 'sigma', 'sigma_err', 'tri', 'dist',
'wrap', 'getdist', 'DBG']
@debug(["distance"])
def __init__(self):
pass
[docs] def setup(self, boxsize, pbc, pos, sigma, tri):
"""
Set up periodic conditions.
Not using __init__ for this allows more decorator flexibility
"""
self.boxsize = boxsize
self.pbc = pbc
self.reset_wrap()
self.posones = ones_like(pos)
self.sigma = sigma
self.sigma_err = self.sigma * 0.7
self.tri = tri
self.dist = empty((pos.shape[-2], pos.shape[-2], 3))
# self.dist_mat = empty((self.nm, self.nm))
return self.distance(pos)
[docs] def reset_wrap(self):
"""Start false."""
if self.pbc:
self.getdist = self._wrapped_dist
self.wrap = self._wrap2
else:
self.getdist = self._dist
self.wrap = self._nowrap
[docs] def setup_full(self, pos):
"""
See setup.
For the full trajectory.
"""
self.getdist = self._dist_full
self.posones = ones_like(pos)
self.dist = empty((*pos.shape[:2], pos.shape[-2], 3))
# self.dist_mat = empty((self.nm, self.nm))
return self.distance_full(pos)
def __call__(self, func_to_decorate):
"""
Periodic Boundary condition decorator.
Position Wrapping
* This will provide and adjust the distances and positions of the atoms
based on boundary condition selection
[REF Computer simulation of liquids Allen and Tildesley pg 30]
"""
@wraps(func_to_decorate)
def wrapper(*args, **kw):
result = func_to_decorate(*args, **kw)
args[0].pos[:] = self.wrap(args[0].pos)
args[0].dist[:], args[0].dist_matrix[:] = self.distance(args[0].pos)
return result
return wrapper
[docs] @debug(['boundaries'])
def distance(self, pos):
"""
Provide a distance matrix for both periodic and non periodic conditions.
Returns
-------
dist: array
All pairs distances n(n+1)/2 by 3
dist_mat: array
distance matrix n by n
only ever used as its inverse so calculated as inverse
"""
self.dist[:] = self.getdist(pos)
# if any(absolute(dist_mat[self.tri]) < self.sigma_err):
# Error.count("W", "Particles overlap")
# print(self.sigma, self.sigma_err, absolute(dist_mat[self.tri]) - self.sigma_err);exit(0)
return self.dist, 1 / maximum(sqrt(c_einsum('ijz,ijz->ij', self.dist, self.dist)), c.EPS_val)
[docs] def distance_full(self, pos):
"""
See distance.
For the full trajectory
"""
self.dist[:] = self.getdist(pos)
return self.dist, 1 / maximum(sqrt(c_einsum('aijz,aijz->aij', self.dist, self.dist)), c.EPS_val)
def _wrap1(self, coords):
self.wrap = self._wrap2
return self.wrapping(coords)
def _wrap2(self, coords):
self.wrap = self._wrap1
return coords
@staticmethod
def _nowrap(coords):
return coords
def _wrapped_dist(self, pos):
"""
Wrap distance every time.
Positions only need to be wrapped once per loop whereas the wrapped distances are always needed
"""
return self.wrapping(self._dist(pos))
def _dist(self, pos):
dummy = self.posones
return c_einsum("iz,jz->ijz", pos, dummy) - c_einsum("iz,jz->jiz", pos, dummy)
def _dist_full(self, pos):
dummy = self.posones
return c_einsum("aiz,ajz->aijz", pos, dummy) - c_einsum("aiz,ajz->ajiz", pos, dummy)
[docs] def wrapping(self, pos):
"""Simplistic wrapping of coordinates."""
return pos - self.boxsize * floor(pos / self.boxsize + 0.5)