Marcello-Sega/pytim

View on GitHub
pytim/surface.py

Summary

Maintainability
C
1 day
Test Coverage
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4

from __future__ import print_function
from abc import ABCMeta, abstractproperty
import numpy as np
from scipy.spatial import Delaunay, cKDTree
from scipy.interpolate import LinearNDInterpolator
from . import utilities
from . import messages
from .observables import LocalReferenceFrame as LocalReferenceFrame
import weakref


class Surface(object):
    """ Everything about the continuum description of surfaces.

        Any implementation must provide the following methods:

        distance()      -> returns an array of relative positions with respect
                           to the interface
        triangulation() -> a triangulated surface, if available
        regular_grid()  -> the surface elevation on a regular grid, if
                           available
        dump()          -> save to disk in format of choice
    """
    __metaclass__ = ABCMeta

    def __init__(self, interface, options=None):
        self.interface = weakref.proxy(interface)
        self.normal = interface.normal
        self.alpha = interface.alpha
        self.options = options
        self.z = self.normal
        try:
            self.xyz = np.roll(np.array([0, 1, 2]), 2 - self.z)
            self.xy = self.xyz[0:2]
        except:
            self.xyz, self.xy = None, None
        try:
            self._layer = options['layer']
        except (TypeError, KeyError):
            self._layer = 0

    @abstractproperty
    def interpolation(self, inp):
        """ returns interpolated position on the surface """
        return utilities.extract_positions(inp)

    @abstractproperty
    def distance(self, inp, *args, **kargs):
        """ returns distance from the surface """
        positions = utilities.extract_positions(inp)
        distance_array = positions[::, 2]
        return distance_array

    @abstractproperty
    def triangulation(self):
        """ return a scipy.spatial.qhull.Delaunay triangulation
            of the surface
        """
        return triangulation

    @abstractproperty
    def regular_grid(self):
        """ returns the points defining the regular grid in 2 dimensions, and
            the elevation values evaluated at the grid points
            (like the input of scipy.interpolate.RegularGridInterpolator )
        """
        return (x, y), elevations

    @abstractproperty
    def dump(self):
        """ save the surface to file """
        return True

    def update_q_vectors(self, box):
        try:
            if np.any(box != self.box):
                self._compute_q_vectors(box)
        except BaseException:
            self._compute_q_vectors(box)

    def _compute_q_vectors(self, box):
        """ Compute the q-vectors compatible with the current box dimensions.

            Calculated quantities:
            q_vectors : two 2D arrays forming the grid of q-values, similar
                        to a meshgrid
            Qxy       : array of the different q-vectors
            Q         : squared module of Qxy with the first element missing
                        (no Q = 0.0)
        """
        self.box = np.roll(box, 2 - self.normal)
        nmax = list(map(int, np.ceil(self.box[0:2] / self.alpha)))
        self.q_vectors = np.mgrid[0:nmax[0], 0:nmax[1]] * 1.0
        self.q_vectors[0] *= 2. * np.pi / box[0]
        self.q_vectors[1] *= 2. * np.pi / box[1]
        self.modes_shape = self.q_vectors[0].shape
        qx = self.q_vectors[0][:, 0]
        qy = self.q_vectors[1][0]
        Qx = np.repeat(qx, len(qy))
        Qy = np.tile(qy, len(qx))
        self.Qxy = np.vstack((Qx, Qy)).T
        self.Q = np.sqrt(np.sum(self.Qxy * self.Qxy, axis=1)[1:])

    @staticmethod
    def _surface_from_modes(points, q_vectors, modes):
        elevation = []
        for point in points:
            dotp = q_vectors[0] * point[0] + q_vectors[1] * point[1]
            phase = np.cos(dotp) + 1.j * np.sin(dotp)
            elevation.append(np.sum((phase * modes).real))
        return np.array(elevation)

    def surface_from_modes(self, points, modes):
        return self._surface_from_modes(points, self.q_vectors, modes)

    def surface_modes(self, points):
        QR = np.dot(self.Qxy, points[:, 0:2].T).T
        # ph[0] are the phases associated to each of the ~ n^2 modes for
        # particle 0.
        # We exclude the zero mode.
        ph = (np.cos(QR) + 1.j * np.sin(QR))[:, 1:]
        z = points[:, 2]
        az = np.mean(z)
        z = z - az
        A = (ph / self.Q)
        z = z
        s = np.dot(np.linalg.pinv(A), z)
        # return the surface modes reshaped into an array
        return np.append(az + 0.j, s / self.Q).reshape(self.modes_shape)

    def triangulate_layer_flat(self, layer=0):
        """Triangulate a layer.

        :param int layer:  (default: 0) triangulate this layer (on both sides\
                           of the interface)
        :return list triangulations:  a list of two Delaunay triangulations,\
                           which are also stored in self.surf_triang
        """
        if layer > len(self.interface._layers[0]):
            raise ValueError(messages.UNDEFINED_LAYER)

        box = self.interface.universe.dimensions[:3][self.xyz]

        upper = (self.interface._layers[0][layer])
        lower = (self.interface._layers[1][layer])
        delta = self.interface.alpha * 4.0 + 1e-6
        # here we rotate the system to have normal along z
        upperpos = utilities.generate_periodic_border(
            upper.positions[:, self.xyz], box, delta, method='2d')[0]
        lowerpos = utilities.generate_periodic_border(
            lower.positions[:, self.xyz], box, delta, method='2d')[0]

        self.surf_triang = [None, None]
        self.trimmed_surf_triangs = [None, None]
        self.triangulation_points = [None, None]
        self.surf_triang[0] = Delaunay(upperpos[:, 0:2])
        self.surf_triang[1] = Delaunay(lowerpos[:, 0:2])
        self.triangulation_points[0] = upperpos[:]
        self.triangulation_points[1] = lowerpos[:]
        self.trimmed_surf_triangs[0] = utilities.trim_triangulated_surface(
            self.surf_triang[0], box)
        self.trimmed_surf_triangs[1] = utilities.trim_triangulated_surface(
            self.surf_triang[1], box)
        return self.surf_triang

    @staticmethod
    def local_env_com(positions, reference_pos, box, nneigh):
        tree = cKDTree(positions, boxsize=box)
        # local_env are the positions of the atoms next to pos
        _, local_env_indices = tree.query(reference_pos, k=nneigh)
        local_env = positions[local_env_indices].copy()
        for k in range(nneigh):
            local_env[:, k, :] = utilities.pbc_compact(local_env[:, k, :],
                                                       reference_pos, box)
        return np.average(local_env, axis=1)

    def _distance_generic(self, inp, symmetry):

        inter = self.interface
        pos = utilities.extract_positions(inp)

        if len(pos) == 0:
            raise ValueError("empty group")
        box = inter.universe.dimensions[:3]
        l1centers = inter.atoms.positions[inter.atoms.layers == 1]
        # tree of the surface triangles' centers
        tree = cKDTree(l1centers, boxsize=box)

        # distances and indices of the surface triangles' centers to pos[]
        dist, ind = tree.query(pos, k=1)

        if not isinstance(inp, np.ndarray):
            # a group has been passed, we know exactly which atoms are
            # surface ones.
            try:
                dist[inp.atoms.layers == 1] = 0.0
                # Warning, the corresponding values of 'ind' will be wrong.
                # If you change this code, check that it won't
                # depend on 'ind'
            except AttributeError:
                raise RuntimeError(
                    "Wrong parameter passed to _distance_generic")
        l1centers_ = l1centers[ind]
        buried = inter.is_buried(pos)
        sign = np.ones(dist.shape[0])
        sign[np.where(buried)[0]] = -1.0
        return sign * dist

        raise ValueError("Incorrect symmetry used for distance calculation")

    def _distance_spherical(self, positions):
        box = self.interface.universe.dimensions[:3]
        cond = self.interface.atoms.layers == 1
        tree = cKDTree(self.interface.atoms.positions[cond], boxsize=box)
        d, i = tree.query(positions, k=1)
        dd, surf_neighs = tree.query(self.interface.atoms.positions[cond][i],
                                     5)
        p = self.interface.atoms.positions[cond]
        center = np.mean(p)

        side = (np.sum((positions - center)**2, axis=1) > np.sum(
            (p - center)**2, axis=1)) * 2 - 1
        return d * side

    def _distance_flat(self, positions):
        box = self.interface.universe.dimensions[:3]

        pos = np.copy(positions)

        cond = np.where(pos[:, 0:2] > box[0:2])
        pos[cond] -= box[cond[1]]
        cond = np.where(pos[:, 0:2] < 0 * box[0:2])
        pos[cond] += box[cond[1]]

        distance = self.interpolation(pos)
        if np.sum(np.isnan(distance)) > 1 + int(0.01 * len(pos)):
            Warning("more than 1% of points have fallen outside"
                    "the convex hull while determining the"
                    "interpolated position on the surface."
                    "Something is wrong.")
        # positive values are outside the surface, negative inside
        cond = np.where(np.isclose(distance, 0., atol=1e-2))
        distance[cond] = 0.0
        return distance

    def _initialize_distance_interpolator_flat(self, layer):
        self._layer = layer
        self.triangulate_layer_flat(layer=self._layer)

        self._interpolator = [None, None]
        for side in [0, 1]:
            self._interpolator[side] = LinearNDInterpolator(
                self.surf_triang[side], self.triangulation_points[side][:, 2])


class SurfaceFlatInterface(Surface):
    def distance(self, inp, *args, **kargs):
        positions = utilities.extract_positions(inp)
        return self._distance_flat(positions)

    def interpolation(self, inp):
        positions = utilities.extract_positions(inp)
        box = self.interface.universe.dimensions[self.z]
        try:
            self.options['from_modes']
            upper_interp = self.surface_from_modes(upper_set, self.modes[0])
            lower_interp = self.surface_from_modes(lower_set, self.modes[1])
        except (TypeError, KeyError):
            self._initialize_distance_interpolator_flat(layer=self._layer)
            upper_interp = self._interpolator[0](positions[:, self.xy])
            lower_interp = self._interpolator[1](positions[:, self.xy])

        d1 = upper_interp - positions[:, self.z]
        d1[np.where(d1 > box / 2.)] -= box
        d1[np.where(d1 < -box / 2.)] += box

        d2 = lower_interp - positions[:, self.z]
        d2[np.where(d2 > box / 2.)] -= box
        d2[np.where(d2 < -box / 2.)] += box

        cond = np.where(np.abs(d1) <= np.abs(d2))[0]
        distance = d2
        distance[cond] = -d1[cond]

        return distance

    def dump(self):
        pass

    def regular_grid(self):
        pass

    def triangulation(self, layer=0):
        return self.triangulate_layer_flat(layer)


class SurfaceGenericInterface(Surface):
    def distance(self, inp, *args, **kargs):
        symmetry = args[0]
        try:
            mode = kargs['mode']
        except:
            mode = 'default'

        if mode == 'default':
            return self._distance_generic(inp, symmetry)

    def interpolation(self, inp):
        pass

    def dump(self):
        pass

    def regular_grid(self):
        pass

    def triangulation(self, layer=0):
        pass


#