# -*- coding: utf-8 -*-
# Copyright (c) 2023, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""Definition of a class for a simulation box.
The simulation box handles the periodic boundaries if needed.
It is typically referenced via the :py:class:`.System` class,
i.e. as ``System.box``.
Important classes defined here
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
BoxBase (:py:class:`.BoxBase`)
Class for a simulation box.
RectangularBox (:py:class:`.RectangularBox`)
Class representing a rectangular simulation box.
TriclinicBox (:py:class:`.TriclinicBox`)
Class representing a triclinic simulation box.
Examples
~~~~~~~~
>>> from pyretis.core.box import create_box
>>> box = create_box(cell=[10, 10, 10], periodic=[True, False, True])
"""
from abc import ABCMeta, abstractmethod
from copy import copy
import logging
import math
import numpy as np
from numpy.linalg import det
from numpy import prod
from pyretis.core.common import compare_objects
logger = logging.getLogger(__name__) # pylint: disable=invalid-name
logger.addHandler(logging.NullHandler())
__all__ = ['create_box', 'check_consistency', 'box_matrix_to_list',
'box_vector_angles', 'angles_from_box_matrix']
[docs]def create_box(low=None, high=None, cell=None, periodic=None):
"""Set up and create a box.
Parameters
----------
low : numpy.array, optional
1D array containing the lower bounds of the cell.
high : numpy.array, optional
1D array containing the higher bounds of the cell.
cell : numpy.array, optional
1D array, a flattened version of the simulation box matrix.
This array is expected to contain 1, 2, 3, 6 or 9 items.
These are the xx, yy, zz, xy, xz, yx, yz, zx, zy elements,
respectively.
periodic : list of boolean, optional
If `periodic[i]` then we should apply periodic boundaries
to dimension `i`.
Returns
-------
out : object like :py:class:`.BoxBase`
The object representing the simulation box.
"""
# If the cell is given, the we should be able to determine
# the length:
if cell is not None:
# make sure the cell does not become an array of objects
cell = np.array(cell, dtype=float)
if len(cell) <= 3:
length = np.array(list(cell))
else: # Use the xx, yy and zz parameters:
length = np.array(list(cell[:3]))
# Determine low, high, length and possible periodic:
low, high, length, periodic = _get_low_high_length(
low, high, length, periodic
)
else: # Here, cell was not given. Try to obtain it from the length:
low, high, length, periodic = _get_low_high_length(
low, high, cell, periodic
)
cell = np.array(list(length))
# We can still have periodic not set:
if periodic is None:
logger.info(
'Periodic settings not given. Assumed True for all directions.'
)
periodic = [True] * len(length)
else:
# If for some reason the periodic settings have wrong length:
if len(periodic) < len(length):
logger.info('Setting missing periodic settings to True.')
for _ in range(len(length) - len(periodic)):
periodic.append(True)
elif len(periodic) > len(length):
logger.error('Too many periodic settings given.')
raise ValueError('Too many periodic settings given.')
# Here, everything should be set:
check_consistency(low, high, length)
# Create the box:
obj = TriclinicBox
if len(cell) <= 3:
obj = RectangularBox
return obj(low, high, length, periodic, cell)
[docs]def check_consistency(low, high, length):
"""Check that given box bounds are consistent.
Parameters
----------
low : numpy.array
The lower bounds for the box.
high : numpy.array
The upper bounds for the box.
length : numpy.array
The lengths of the box.
"""
length_re = high - low
if any(i <= 0 for i in length_re):
logger.error('Check box settings. Found high <= low.')
raise ValueError('Incorrect box: high <= low.')
if not all(np.isclose(length, length_re)):
logger.error('Check box settings: length != high - low.')
raise ValueError('Incorrect box: length != high - low.')
[docs]def _get_low_high_length(low, high, length, periodic):
"""Determine box cell parameters from input.
This method will consider the following cases:
1) We are given low, high and length.
2) We are given high and length, and determine the low values.
3) We are given low and length, and determine the high values.
4) We are given length, assume low to be zero and determine high.
5) We are given low and high, and determine the length.
6) We are given just high, assume low to be zero and
determine the length.
7) We are just given low, and assume high and length to be infinite.
8) We are given none of the values, and assume all to infinite.
Parameters
----------
low : numpy.array or None
The lower bounds for the box.
high : numpy.array or None
The upper bounds for the box.
length : numpy.array or None
The lengths of the box.
periodic : list of boolean or None
We will assume a periodic box for dimensions where
this list is True.
Returns
-------
out[0] : numpy.array
The updated lower bounds for the box.
out[1] : numpy.array
The updated upper bounds for the box.
out[2] : numpy.array
The updated lengths of the box.
out[3] : list of boolean
The updated periodic settings for the box.
"""
case = (length is not None, low is not None, high is not None)
if low is not None:
low = np.array(low)
if high is not None:
high = np.array(high)
if length is not None:
length = np.array(length)
if case == (True, True, True):
# 1) We have given length, low and high.
pass
elif case == (True, False, True):
# 2) Length & high has been given, just determine low:
low = high - length
elif case == (True, True, False):
# 3) Length and low was given, determine high:
high = low + length
elif case == (True, False, False):
# 4) Length is given, set low to 0 and high to low + length:
low = np.zeros_like(length)
high = low + length
elif case == (False, True, True):
# 5) Low and high is given, determine length:
length = high - low
elif case == (False, False, True):
# 6) High is given, assume low and determine length:
low = np.zeros_like(high)
length = high - low
elif case == (False, True, False):
# 7) Low given. High and length to be determined.
# This is not enough info, so we assume an infinite box:
length = float('inf') * np.ones_like(low)
high = float('inf') * np.ones_like(low)
elif case == (False, False, False):
# Not much info is given. We let the box be similar
# in shape to the input periodic settings:
if periodic is None:
logger.info(
'Too few settings for the box is given. A 1D box is assumed.'
)
periodic = [False]
low = np.array([-float('inf') for _ in periodic])
high = float('inf') * np.ones_like(low)
length = float('inf') * np.ones_like(low)
return low, high, length, periodic
def array_to_box_matrix(cell):
"""Return a box matrix corresponding to a cell array.
Parameters
----------
cell : list or numpy.array
An (1D) array containing 1, 2, 3, 6 or 9 items. These are
the xx, yy, zz, xy, xz, yx, yz, zx, zy elements. Setting
x = 0, y = 1 and z = 2 will give the indices in the matrix,
e.g. yx -> (1, 0) will correspond to the item in row 1 and
column 0.
Returns
-------
box : numpy.array (2D)
The box vector on matrix form.
"""
if len(cell) == 1:
return 1.0 * np.array([cell[0]])
if len(cell) == 2:
return 1.0 * np.array([[cell[0], 0.0],
[0.0, cell[1]]])
if len(cell) == 3:
return 1.0 * np.array([[cell[0], 0.0, 0.0],
[0.0, cell[1], 0.0],
[0.0, 0.0, cell[2]]])
if len(cell) == 6:
return 1.0 * np.array([[cell[0], cell[3], cell[4]],
[0.0, cell[1], cell[5]],
[0.0, 0.0, cell[2]]])
if len(cell) == 9:
return 1.0 * np.array([[cell[0], cell[3], cell[4]],
[cell[5], cell[1], cell[6]],
[cell[7], cell[8], cell[2]]])
logger.error(
'%d box parameters given, need 1, 2, 3, 6, or 9.', len(cell)
)
raise ValueError('Incorrect number of box-parameters!')
[docs]def box_matrix_to_list(matrix, full=False):
"""Return a list representation of the box matrix.
This method ensures correct ordering of the elements for PyRETIS:
``xx, yy, zz, xy, xz, yx, yz, zx, zy``.
Parameters
----------
matrix : numpy.array
A matrix (2D) representing the box.
full : boolean, optional
Return a full set of parameters (9) if set to True. If False,
and we need 3 or fewer parameters (i.e. the other 6 are zero)
we will only return the 3 non-zero ones.
Returns
-------
out : list
A list with the box-parametres.
"""
if matrix is None:
return None
if np.count_nonzero(matrix) <= 3 and not full:
return [matrix[0, 0], matrix[1, 1], matrix[2, 2]]
return [matrix[0, 0], matrix[1, 1], matrix[2, 2],
matrix[0, 1], matrix[0, 2], matrix[1, 0],
matrix[1, 2], matrix[2, 0], matrix[2, 1]]
[docs]def _cos(angle):
"""Return cosine of an angle.
We also check if the angle is close to 90.0 and if so, we return
a zero.
Parameters
----------
angle : float
The angle in degrees.
Returns
-------
out : float
The cosine of the angle.
"""
if math.isclose(angle, 90.):
return 0.
return math.cos(math.radians(angle))
[docs]def box_vector_angles(length, alpha, beta, gamma):
"""Obtain the box matrix from given lengths and angles.
Parameters
----------
length : numpy.array
1D array, the box-lengths on form ``[a, b, c]``.
alpha : float
The alpha angle, in degrees.
beta : float
The beta angle, in degrees.
gamma : float
The gamma angle, in degrees.
Returns
-------
out : numpy.array, 2D
The (upper triangular) box matrix.
"""
box_matrix = np.zeros((3, 3))
cos_alpha = _cos(alpha)
cos_beta = _cos(beta)
cos_gamma = _cos(gamma)
box_matrix[0, 0] = length[0]
box_matrix[0, 1] = length[1] * cos_gamma
box_matrix[0, 2] = length[2] * cos_beta
box_matrix[1, 1] = math.sqrt(length[1]**2 - box_matrix[0, 1]**2)
box_matrix[1, 2] = (length[1] * length[2] * cos_alpha -
box_matrix[0, 1] * box_matrix[0, 2]) / box_matrix[1, 1]
box_matrix[2, 2] = math.sqrt(length[2]**2 - box_matrix[0, 2]**2 -
box_matrix[1, 2]**2)
return box_matrix
[docs]def angles_from_box_matrix(box_matrix):
"""Obtain angles and lengths from a given box matrix.
Parameters
----------
box_matrix : np.array
The box matrix (2D).
Returns
-------
length : numpy.array
1D array, the box-lengths on form ``[a, b, c]``.
alpha : float
The alpha angle, in degrees.
beta : float
The beta angle, in degrees.
gamma : float
The gamma angle, in degrees.
"""
length = [
box_matrix[0, 0],
math.sqrt(box_matrix[1, 1]**2 + box_matrix[0, 1]**2),
math.sqrt(box_matrix[2, 2]**2 + box_matrix[0, 2]**2 +
box_matrix[1, 2]**2),
]
cos_alpha = (
(box_matrix[0, 1] * box_matrix[0, 2] +
box_matrix[1, 1] * box_matrix[1, 2]) / (length[1] * length[2])
)
cos_beta = box_matrix[0, 2] / length[2]
cos_gamma = box_matrix[0, 1] / length[1]
alpha = math.degrees(math.acos(cos_alpha))
beta = math.degrees(math.acos(cos_beta))
gamma = math.degrees(math.acos(cos_gamma))
return np.array(length), alpha, beta, gamma
class BoxBase(metaclass=ABCMeta):
"""Class for a generic simulation box.
This class defines a generic simulation box.
Attributes
----------
low : numpy.array
1D array containing the lower bounds of the cell.
high : numpy.array
1D array containing the higher bounds of the cell.
length : numpy.array
1D array containing the length of the sides of the
simulation box.
ilength : numpy.array
1D array containing the inverse box lengths for the
simulation box.
periodic : list of boolean
If `periodic[i]` then we should apply periodic boundaries
to dimension `i`.
box_matrix : numpy.array
2D matrix, representing the simulation cell.
cell : numpy.array
1D array representing the simulation cell (flattened
version of the 2D box matrix).
"""
def __init__(self, low, high, length, periodic, cell):
"""Initialise the BoxBase class."""
self.low = low
self.high = high
self.length = length
self.periodic = periodic
self.cell = cell
# Create box matrix from the given cell:
self.box_matrix = array_to_box_matrix(self.cell)
self.ilength = 1.0 / self.length
self.dim = len(self.length)
def update_size(self, new_size):
"""Update the box size.
Parameters
----------
new_size : list, tuple, numpy.array, or other iterable
The new box size.
"""
if new_size is None:
logger.warning(
'Box update ignored: Tried to update with empty size!'
)
else:
try:
size = new_size.size
except AttributeError:
size = len(new_size)
if size <= 3:
if size == len(self.cell):
for i in range(self.dim):
self.length[i] = new_size[i]
self.high[i] = self.low[i] + new_size[i]
self.cell[i] = new_size[i]
self.ilength = 1.0 / self.length
else:
try:
self.box_matrix = array_to_box_matrix(new_size)
self.cell = list(new_size)
self.length = np.array([float(i) for i in self.cell[:3]])
self.high = self.low + self.length
self.ilength = 1.0 / self.length
except ValueError:
logger.critical('Box update failed!')
def bounds(self):
"""Return the boundaries of the box (low, high) as an array."""
return list(zip(self.low, self.high))
@abstractmethod
def calculate_volume(self):
"""Return the volume of the box."""
return
@abstractmethod
def pbc_coordinate_dim(self, pos, dim):
"""Apply periodic boundaries to a selected dimension only.
For the given positions, this function will apply periodic
boundary conditions to one dimension only. This can be useful
for instance in connection with order parameters.
Parameters
----------
pos : float
Coordinate to wrap.
dim : int
This selects the dimension to consider.
"""
return
@abstractmethod
def pbc_wrap(self, pos):
"""Apply periodic boundaries to the given position.
Parameters
----------
pos : nump.array
Positions to apply periodic boundaries to.
Returns
-------
out : numpy.array, same shape as parameter `pos`
The periodic-boundary wrapped positions.
"""
return
@abstractmethod
def pbc_dist_matrix(self, distance):
"""Apply periodic boundaries to a distance matrix/vector.
Parameters
----------
distance : numpy.array
The distance vectors.
Returns
-------
out : numpy.array, same shape as parameter `distance`
The periodic-boundary-wrapped distances.
"""
return
@abstractmethod
def pbc_dist_coordinate(self, distance):
"""Apply periodic boundaries to a distance.
This will apply periodic boundaries to a distance. Note that the
distance can be a vector, but not a matrix of distance vectors.
Parameters
----------
distance : numpy.array with shape `(self.dim,)`
A distance vector.
Returns
-------
out : numpy.array, same shape as parameter `distance`
The periodic-boundary wrapped distance vector.
"""
return
def print_length(self, fmt=None):
"""Return a string with box lengths. Can be used for output."""
if fmt is None:
return ' '.join((f'{i}' for i in self.cell))
return ' '.join((fmt.format(i) for i in self.cell))
def restart_info(self):
"""Return a dictionary with restart information."""
info = {
'length': self.length,
'periodic': self.periodic,
'low': self.low,
'high': self.high,
'cell': self.cell,
}
return info
def copy(self):
"""Return a copy of the box.
Returns
-------
out : object like :py:class:`.BoxBase`
A copy of the box.
"""
box_copy = self.__class__(
np.copy(self.low),
np.copy(self.high),
np.copy(self.length),
copy(self.periodic),
np.copy(self.cell)
)
return box_copy
def load_restart_info(self, info):
"""Read the restart information."""
self.length = info.get('length')
self.periodic = info.get('periodic')
self.low = info.get('low')
self.high = info.get('high')
self.cell = info.get('cell')
def __str__(self):
"""Return a string describing the box.
Returns
-------
out : string
String with the type of box, the extent of the box and
information about the periodicity.
"""
boxstr = []
if len(self.cell) <= 3:
boxstr.append('Orthogonal box:')
else:
boxstr.append('Triclinic box:')
for i, periodic in enumerate(self.periodic):
low = self.low[i]
high = self.high[i]
msg = 'Dim: {}, Low: {}, high: {}, periodic: {}'
boxstr.append(msg.format(i, low, high, periodic))
cell = self.print_length()
boxstr.append(f'Cell: {cell}')
return '\n'.join(boxstr)
def __eq__(self, other):
"""Compare two box objects."""
attrs = {'low', 'high', 'length', 'ilength', 'box_matrix', 'cell',
'periodic', 'dim'}
numpy_attrs = {'low', 'high', 'length', 'ilength', 'box_matrix',
'cell'}
return compare_objects(self, other, attrs, numpy_attrs)
def __ne__(self, other):
"""Compare two box objects."""
return not self == other
class RectangularBox(BoxBase):
"""An orthogonal box."""
def calculate_volume(self):
"""Calculate the volume of the box.
Returns
-------
out : float
The volume of the box.
"""
return prod(self.length)
def pbc_coordinate_dim(self, pos, dim):
"""Apply periodic boundaries to a selected dimension only.
For the given positions, this function will apply periodic
boundary conditions to one dimension only. This can be useful
for instance in connection with order parameters.
Parameters
----------
pos : float
Coordinate to wrap around.
dim : int
This selects the dimension to consider.
"""
if self.periodic[dim]:
low, length = self.low[dim], self.length[dim]
ilength = self.ilength[dim]
relpos = pos - low
delta = relpos
if relpos < 0.0 or relpos >= length:
delta = relpos - np.floor(relpos * ilength) * length
return delta + low
return pos
def pbc_wrap(self, pos):
"""Apply periodic boundaries to the given position.
Parameters
----------
pos : nump.array
Positions to apply periodic boundaries to.
Returns
-------
out : numpy.array, same shape as parameter `pos`
The periodic-boundary wrapped positions.
"""
pbcpos = np.zeros(pos.shape)
for i, periodic in enumerate(self.periodic):
if periodic:
low = self.low[i]
length = self.length[i]
ilength = self.ilength[i]
relpos = pos[:, i] - low
delta = np.where(
np.logical_or(relpos < 0.0, relpos >= length),
relpos - np.floor(relpos * ilength) * length,
relpos
)
pbcpos[:, i] = delta + low
else:
pbcpos[:, i] = pos[:, i]
return pbcpos
def pbc_dist_matrix(self, distance):
"""Apply periodic boundaries to a distance matrix/vector.
Parameters
----------
distance : numpy.array
The distance vectors.
Returns
-------
out : numpy.array, same shape as the `distance` parameter
The pbc-wrapped distances.
Note
----
This will modify the given input matrix inplace. This can be
changed by setting ``pbcdist = np.copy(distance)``.
"""
pbcdist = distance
for i, (periodic, length, ilength) in enumerate(zip(self.periodic,
self.length,
self.ilength)):
if periodic:
dist = pbcdist[:, i]
high = 0.5 * length
k = np.where(np.abs(dist) >= high)[0]
dist[k] -= np.rint(dist[k] * ilength) * length
return pbcdist
def pbc_dist_coordinate(self, distance):
"""Apply periodic boundaries to a distance.
This will apply periodic boundaries to a distance. Note that the
distance can be a vector, but not a matrix of distance vectors.
Parameters
----------
distance : numpy.array with shape `(self.dim,)`
A distance vector.
Returns
-------
out : numpy.array, same shape as the `distance` parameter
The periodic-boundary wrapped distance vector.
"""
pbcdist = np.zeros(distance.shape)
for i, (periodic, length, ilength) in enumerate(zip(self.periodic,
self.length,
self.ilength)):
if periodic and np.abs(distance[i]) > 0.5*length:
pbcdist[i] = (distance[i] -
np.rint(distance[i] * ilength) * length)
else:
pbcdist[i] = distance[i]
return pbcdist
class TriclinicBox(BoxBase):
"""This class represents a triclinic box."""
def calculate_volume(self):
"""Calculate and return the volume of the box.
Returns
-------
out : float
The volume of the box.
"""
return det(self.box_matrix)
def pbc_coordinate_dim(self, pos, dim):
"""Apply periodic boundaries to a selected dimension only."""
raise NotImplementedError
def pbc_wrap(self, pos):
"""Apply periodic boundaries to the given position."""
raise NotImplementedError
def pbc_dist_matrix(self, distance):
"""Apply periodic boundaries to a distance matrix/vector."""
raise NotImplementedError
def pbc_dist_coordinate(self, distance):
"""Apply periodic boundaries to a distance."""
raise NotImplementedError
def box_from_restart(restart):
"""Create a box from restart settings.
Parameters
----------
restart : dict
A dictionary with restart settings.
Returns
-------
box : object like :py:class:`.BoxBase`
The box created from the restart settings.
"""
restart_box = restart.get('box', None)
if restart_box is None:
logger.info('No box created from restart settings.')
return None
box = create_box(
low=restart_box.get('low'),
high=restart_box.get('high'),
cell=restart_box.get('cell'),
periodic=restart_box.get('periodic')
)
return box