# -*- coding: utf-8 -*-
import numpy
from copy import deepcopy
from kvazar.core.lib import libmd
try:
from mpi4py import MPI
comm = MPI.comm
except:
comm = None
class BoxcellException(Exception):
def __init__(self, message):
super(BoxcellException, self).__init__(message)
[docs]class BoxCell:
"""
It's a class that realizes periodic boundary conditions.
"""
def __init__(self, coords, r_in, r_out, box):
self.r_in = float(r_in)
self.r_out = float(r_out)
self.delta = r_out - r_in
self.n_cells = numpy.array(map(int, box[:, 2] / r_out), numpy.int32)
self.r_cells = numpy.array(box[:,2] / self.n_cells, numpy.float)
self.box = box
self.update_prev_coords(coords)
[docs] def update_prev_coords(self, coords):
"""
Rewriting previous coordinates.
"""
self.prev_coords = deepcopy(coords)
[docs] def time_to_update(self, coords, a, b):
"""
If coordinates have changed significantly then cutoff lists should be
updated.
"""
if comm:
max_len_ab = libmd.max_pbc_dist(self.prev_coords, cur_coords, a, b)
max_len = comm.allreduce(max_len_ab, MPI.MAX)
return self.delta < 2 * max_len
else:
return self.delta < 2 * libmd.max_pbc_dist(self.prev_coords, coords, self.box, 0, coords.shape[0])
[docs] def build_cells(self, coords):
"""
This method checks sizes of structure and raises exception if structure
is bigger than periodic box. Else it returns list of atom names for each
part of periodic box.
"""
length = numpy.prod(self.n_cells)
n_cells = self.n_cells
r_cells = self.r_cells
cells = [[] for _ in xrange(length)]
xlim = coords[:,0].min(), coords[:,0].max()
ylim = coords[:,1].min(), coords[:,1].max()
zlim = coords[:,2].min(), coords[:,2].max()
checkx = self.box[0, 0] < xlim[0] and self.box[0, 1] > xlim[1]
checky = self.box[1, 0] < ylim[0] and self.box[1, 1] > ylim[1]
checkz = self.box[2, 0] < zlim[0] and self.box[2, 1] > zlim[1]
if not (checkx and checky and checkz):
message = "Too large molecular structure for this periodic box.\n"
message += "BOX_X : [%s, %s]\n" % (self.box[0, 0], self.box[0, 1])
message += "BOX_Y : [%s, %s]\n" % (self.box[1, 0], self.box[1, 1])
message += "BOX_Z : [%s, %s]\n" % (self.box[2, 0], self.box[2, 1])
message += "STRUCTURE_X : [%s, %s]\n" % xlim
message += "STRUCTURE_Y : [%s, %s]\n" % ylim
message += "STRUCTURE_Z : [%s, %s]\n" % zlim
raise BoxcellException(message)
ijk = numpy.int32((coords - self.box[:,0]) / r_cells)
mulind = numpy.array([n_cells[1] * n_cells[2], n_cells[2], 1], numpy.int32)
ijk *= mulind
ijk = ijk.sum(axis = 1)
[cells[j].append(i) for i, j in enumerate(ijk)]
return cells
def neighbours_for_crd(self, coords, crd, cells = None):
if not cells:
cells = self.build_cells(coords)
ijk = numpy.int32(crd - self.box[:,0] / self.r_cells)
index = ijk.sum(axis = 1)
neigh = [[] for _ in xrange(len(index))]
libmd.make_cutoff_nb_lists_twostruct(neigh, cells, crd, coords, self.r_cells, self.n_cells, self.box, self.r_out, 0, len(index))
return neigh
[docs] def calc_neighbours(self, coords, a, b):
"""
It builds periodic box and then update cutoff lists, if needed.
"""
cells = self.build_cells(coords)
neighbours = [[] for _ in range(coords.shape[0])]
libmd.make_cutoff_nb_lists(neighbours, cells, coords, self.r_cells, self.n_cells, self.box, self.r_out, a, b)
return neighbours