Merge pull request #20 from LaurentRDC/develop

Develop

docs/source/api.rst

@@ -38,13 +38,17 @@ Affine Transforms
 
 Structure
 =========
-.. autoclass:: skued.structure.Atom
+Handling crystal structure information is crucial for many data analysis and modelling tasks.
+See the :ref:`Structure tutorial <structure_tutorial>` for some examples on how to use the following
+classes.
+
+.. autoclass:: skued.Atom
    :members:
 
-.. autoclass:: skued.structure.Lattice
+.. autoclass:: skued.Lattice
    :members:
 
-.. autoclass:: skued.structure.Crystal
+.. autoclass:: skued.Crystal
    :members:
 
 Simulation

docs/source/tutorials/structure.rst

@@ -24,7 +24,7 @@ according to any affine transform.
 
 To create an atom, simply provide its element and coordinates::
 	
-	from skued.structure import Atom
+	from skued import Atom
 
 	copper = Atom(element = 'Cu', coords = [0,0,0])
 
@@ -80,7 +80,7 @@ Constructing a :class:`Crystal` object
 --------------------------------------
 Creating a :class:`Crystal` object can be done most easily from a Crystal Information File (CIF, .cif)::
 	
-	from skued.structure import Crystal
+	from skued import Crystal
 
 	TiSe2 = Crystal.from_cif('tise2.cif')
 
@@ -93,7 +93,7 @@ constructed like so::
 Protein DataBank files are even easier to handle; simply provide the 4-letter identification code
 and the structure file will be taken care of by scikit-ued::
 	
-	bacteriorhodopsin = Crystal.from_pdb('1fbb')
+	hemoglobin = Crystal.from_pdb('1gxz')
 
 Another convenient way to construct a :class:`Crystal` is through the `Crystallography Open Database <http://www.crystallography.net/cod/>`_::
 
@@ -119,7 +119,7 @@ not need to be provided. The symmetry operators must be expressed in the reduced
 As an example, let's create the simplest crystal structure known: 
 `alpha-Polonium (simple cubic) <https://en.wikipedia.org/wiki/Polonium#Solid_state_form>`_::
 	
-	from skued.structure import Crystal, Atom
+	from skued import Crystal, Atom
 	import numpy as np
 
 	lattice_vectors = 3.35 * np.eye(3)
@@ -141,12 +141,6 @@ The :class:`Crystal` object provides some interfaces for easy structure manipula
 Note that iterating over the :attr:`crystal.atoms` attribute may or may not be equivalent to 
 :data:`iter(crystal)`, due to the way symmetry operators are defined.
 
-:class:`Crystal` objects also provide interoperability with :mod:`spglib`::
-
-	import spglib
-
-	spglib.get_symmetry_dataset(cell = graphite.spglib_cell, symprec = 1e-4)
-
 Lattice vectors and reciprocal space
 -------------------------------------
 Once a :class:`Crystal` object is ready, you can manipulate the lattice parameters via the underlying :class:`Lattice`
@@ -166,6 +160,39 @@ The unit cell volume (and by extensions, density) is also accessible:
 	vol = graphite.volume
 	density = vol/len(graphite)
 
+Space-group Information
+-----------------------
+Thanks to `spglib <http://atztogo.github.io/spglib/>`_, we can get symmetry and space-group information 
+from a :class:`Crystal` instance::
+
+	from skued import Crystal
+	
+	gold = Crystal.from_database('Au')
+	spg_info = gold.spacegroup_info()
+
+In the above example, :data:`spg_info` is a dictionary with the following four keys:
+
+* ``'international_symbol'``: International Tables of Crystallography space-group symbol (short);
+
+* ``'international_full'``: International Tables of Crystallography space-group full symbol;
+
+* ``'hall_symbol'`` : Hall symbol;
+
+* ``'pointgroup'`` : International Tables of Crystallography point-group;
+
+* ``'international_number'`` : International Tables of Crystallography space-group number (between 1 and 230);
+
+* ``'hall_number'`` : Hall number (between 1 and 531).
+
+You can get even more information by using :mod:`spglib` functions directly::
+
+	from spglib import get_symmetry_dataset
+
+	all_the_info = get_symmetry_dataset(gold.spglib_cell)
+
+The content of the :data:`all_the_info` dictionary is documented `here <http://atztogo.github.io/spglib/python-spglib.html#get-symmetry-dataset>`_.
+Many of :mod:`spglib`'s routines can be used with :attr:`Crystal.spglib_cell`.
+
 Scattering utilities
 --------------------
 :class:`Crystal` objects have a few methods that make life easier when dealing with scattering data and modeling.
@@ -195,42 +222,4 @@ Static structure factor calculation is also possible, both for a single reflecti
 	SF = graphite.structure_factor_miller(h, k, l)
 	SF.shape == h.shape 	# True
 
-Atomic potential
-----------------
-The crystal electrostatic potential (the scattering potential leading to electron diffraction) can be
-computed from a :class:`Crystal`::
-
-	import numpy as np
-	import matplotlib.pyplot as plt
-	from skued.structure import graphite
-
-	xx, yy = np.meshgrid(np.linspace(-1, 1, num = 100), 
-	                     np.linspace(-1, 1, num = 100))
-	zz = np.zeros_like(xx)
-
-	potential = graphite.potential(xx, yy, zz)
-	plt.imshow(es_potential)
-
-After plot formatting:
-
-.. plot::
-
-	import numpy as np
-	import matplotlib.pyplot as plt
-	from skued.structure import graphite
-	xx, yy = np.meshgrid(np.linspace(-5, 5, num = 512), 
-						 np.linspace(-5, 5, num = 512))
-	zz = np.zeros_like(xx)
-	potential = graphite.potential(xx, yy, zz)
-	plt.title('Electrostatic potential of graphite (log-scale)')
-	plt.imshow(np.log(1 + 0.01*potential), extent = [xx.min(), xx.max(), yy.min(), yy.max()])
-	plt.ylabel('x-direction ($\AA$)')
-	plt.xlabel('y-direction ($\AA$)')
-	plt.show()
-
-Note that while the :data:`graphite` crystal only has four atoms in its unitcell, the :meth:`graphite.potential` method
-will process the meshes :data:`xx`, :data:`yy`, and :data:`zz` through the :func:`skued.minimum_image_distance` function
-to implement the **minimum-image distance convention** for periodic boundary conditions All this to say that the potential 
-is computed for a seemingly-infinite crystal.
-
 :ref:`Return to Top <structure_tutorial>`

skued/__init__.py

@@ -2,7 +2,7 @@
 __author__ = 'Laurent P. René de Cotret'
 __email__ = 'laurent.renedecotret@mail.mcgill.ca'
 __license__ = 'MIT'
-__version__ = '0.4.7' # TODO: automatic versioning?
+__version__ = '0.4.8' # TODO: automatic versioning?
 
 from .affine import (affine_map, change_basis_mesh, change_of_basis, is_basis,
                      is_rotation_matrix, minimum_image_distance,
@@ -13,4 +13,5 @@
 from .parallel import pmap, preduce
 from .plot_utils import spectrum_colors, rgb_sweep
 from .quantities import electron_wavelength, interaction_parameter, lorentz
+from .structure import Crystal, Atom, Lattice
 from .voigt import gaussian, lorentzian, pseudo_voigt

skued/baseline/algorithms.py

@@ -253,17 +253,12 @@ def baseline_dt(array, max_iter, level = None, first_stage = 'sym6', wavelet = '
 	background_regions : iterable, optional
 		Indices of the array values that are known to be purely background. Depending
 		on the dimensions of array, the format is different:
-		
-		``array.ndim == 1``
-			`background_regions` is a list of ints (indices) or slices
-			E.g. >>> background_regions = [0, 7, 122, slice(534, 1000)]
-			
-		``array.ndim == 2``
-			`background_regions` is a list of tuples of ints (indices) or tuples of slices
-			E.g. >>> background_regions = [(14, 19), (42, 99), (slice(59, 82), slice(81,23))]
-			
-		Default is empty list.
 
+		* 1D signal: `background_regions` is a list of ints (indices) or slices, e.g. ``[0, 7, slice(534, 1000)]``.
+
+		* 2D signal: `background_regions` is a list of tuples of ints (indices) or tuples of slices, e.g. ``[(14, 19), (slice(59, 82), slice(81,23))]``.
+
+		Default is empty list.
 	mask : `~numpy.ndarray`, dtype bool, optional
 		Mask array that evaluates to True for pixels that are invalid. 
 	mode : str, optional
@@ -304,20 +299,15 @@ def baseline_dwt(array, max_iter, level = None, wavelet = 'sym6', background_reg
 	wavelet : PyWavelet.Wavelet object or str, optional
 		Wavelet with which to perform the algorithm. See PyWavelet documentation
 		for available values. Default is 'sym6'.
-	background_regions : list, optional
+	background_regions : iterable, optional
 		Indices of the array values that are known to be purely background. Depending
 		on the dimensions of array, the format is different:
-		
-		``array.ndim == 1``
-			`background_regions` is a list of ints (indices) or slices
-			E.g. >>> background_regions = [0, 7, 122, slice(534, 1000)]
-			
-		``array.ndim == 2``
-			`background_regions` is a list of tuples of ints (indices) or tuples of slices
-			E.g. >>> background_regions = [(14, 19), (42, 99), (slice(59, 82), slice(81,23))]
-			
-		Default is empty list.
 
+		* 1D signal: `background_regions` is a list of ints (indices) or slices, e.g. ``[0, 7, slice(534, 1000)]``.
+
+		* 2D signal: `background_regions` is a list of tuples of ints (indices) or tuples of slices, e.g. ``[(14, 19), (slice(59, 82), slice(81,23))]``.
+
+		Default is empty list.
 	mask : `~numpy.ndarray`, dtype bool, optional
 		Mask array that evaluates to True for pixels that are invalid. Useful to determine which pixels are masked
 		by a beam block.

skued/simulation/__init__.py

@@ -2,4 +2,5 @@
 """
 Diffraction simulation
 """
-from .powdersim import powdersim
+from .powdersim import powdersim
+from .potential import electrostatic, pelectrostatic

skued/simulation/potential.py

@@ -0,0 +1,104 @@
+# -*- coding: utf-8 -*-
+"""
+Electrostatic potential simulation
+==================================
+"""
+from functools import partial
+from .. import chunked, minimum_image_distance, repeated_array
+from scipy.special import k0 as bessel
+import numpy as np
+from numpy import pi
+
+m = 9.109*10**(-31)     #in kg
+a0 = 0.5291             #in Angs
+e = 14.4                #Volt*Angstrom
+
+def electrostatic(crystal, x, y, z):
+    """
+    Electrostatic potential from a crystal calculated on a real-space mesh, 
+    assuming an infinite crystal.
+
+    Parameters
+    ----------
+    crystal : skued.Crystal
+        
+    x, y, z : `~numpy.ndarray`
+        Real space coordinates mesh. 
+    
+    Returns
+    -------
+    potential : `~numpy.ndarray`, dtype float
+        Linear superposition of atomic potential [V*Angs]
+    
+    See also
+    --------
+    pelectrostatic
+        Projected electrostatic potential from a crystal
+    """
+    # TODO: split xx and yy into smalled non-repeating unit
+    # TODO: multicore
+    # TODO: pre-load scattering params _a, _b, _c, and _d into an appropriate shape
+
+    potential = np.zeros_like(x, dtype = np.float)
+    r = np.zeros_like(x, dtype = np.float)
+    for atom in crystal:
+        ax, ay, az = atom.xyz(crystal)
+        r[:] = minimum_image_distance(x - ax, y - ay, z - az, 
+                                        lattice = crystal.lattice_vectors)
+        potential += atom.potential(r)
+
+    # Due to sampling, x,y, and z might pass through the center of atoms
+    # Replace np.inf by the next largest value
+    m = potential[np.isfinite(potential)].max()
+    potential[np.isinf(potential)] = m
+    return potential
+
+def pelectrostatic(crystal, xx, yy, bounds = None):
+    """
+    Projected electrostatic potential from a crystal calculated on a real-space mesh, 
+    assuming an infinite crystal. Projection axis is defined as the z-axis. To project the potential
+    along a different axis, the crystal can be rotated with ``Crystal.transform``.
+
+    Parameters
+    ----------
+    crystal : skued.Crystal
+        
+    xx, yy:  `~numpy.ndarray`
+        Real space coordinates mesh. 
+    
+    Returns
+    -------
+    potential : `~numpy.ndarray`, dtype float
+        Linear superposition of atomic potential [V*Angs]
+    """
+    # TODO: split xx and yy into smalled non-repeating unit
+    # TODO: multicore
+    # TODO: pre-load scattering params _a, _b, _c, and _d into an appropriate shape
+
+    if bounds:
+        min_z, max_z = min(bounds), max(bounds)
+        atoms = (atom for atom in iter(crystal) if min_z <= atom.xyz(crystal)[2] < max_z)
+    else:
+        atoms = iter(crystal)
+
+    potential = np.zeros_like(xx)
+    zz = np.zeros_like(xx)
+    for atom in atoms:
+        xa, ya, _ = atom.xyz(crystal)
+        r = minimum_image_distance(xx - xa, yy - ya, zz, lattice = np.array(crystal.lattice_vectors)).reshape(-1,1)
+        potential += np.sum( 2*atom._a*bessel(2*pi*r*np.sqrt(atom._b)) + (atom._c/atom._d) * np.exp( -(r*pi)**2 / atom._d), axis = -1).reshape(xx.shape)
+    potential *= 2 * a0 * e * (pi**2)
+
+    # Due to sampling, x,y, and z might pass through the center of atoms
+    # Replace n.inf by the next largest value
+    potential[np.isinf(potential)] = np.nan
+    potential[np.isnan(potential)] = np.nanmax(potential)
+    return potential
+
+def _proj_elec_atm(atom, xx, yy, lattice):
+    potential = np.zeros_like(xx, dtype = np.float)[:,:,None]
+    xa, ya, _ = tuple(atom.xyz(lattice))
+    r = np.zeros_like(xx)[:,:,None, None]
+    r[:,:, 0, 0] = minimum_image_distance(xx - xa, yy - ya, np.zeros_like(xx), lattice)
+    potential = np.sum( 2*atom._a*bessel(2*pi*r*np.sqrt(atom._b)) + (atom._c/atom._d) * np.exp( -(r*pi)**2 / atom._d), axis = -1)
+    return 2*a0*e*(pi**2)*np.squeeze(potential)   

skued/simulation/tests/test_potential.py

@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+from .. import electrostatic, pelectrostatic
+from ...structure import graphite
+from copy import deepcopy
+import numpy as np
+import unittest
+
+class TestElectrostatic(unittest.TestCase):
+
+	def setUp(self):
+		self.crystal = deepcopy(graphite)
+
+	def test_return_shape(self):
+		""" Test that the return shape of pelectrostatic is the same as input arrays """
+		xx, yy, zz = np.meshgrid(np.linspace(-10, 10, 16), np.linspace(-10, 10, 16), np.linspace(-10, 10, 16))
+		potential = electrostatic(graphite, xx, yy, zz)
+
+		self.assertSequenceEqual(xx.shape, potential.shape)
+
+	def test_side_effects(self):
+		""" Test that mesh arrays are not written to in pelectrostatic """
+		xx, yy, zz = np.meshgrid(np.linspace(-10, 10, 16), np.linspace(-10, 10, 16), np.linspace(-10, 10, 16))
+
+		xx.setflags(write = False)
+		yy.setflags(write = False)
+		zz.setflags(write = False)
+
+		potential = electrostatic(graphite, xx, yy, zz)
+
+class TestPElectrostatic(unittest.TestCase):
+
+	def setUp(self):
+		self.crystal = deepcopy(graphite)
+
+	def test_return_shape(self):
+		""" Test that the return shape of pelectrostatic is the same as input arrays """
+		xx, yy = np.meshgrid(np.linspace(-10, 10, 32), np.linspace(-10, 10, 32))
+		potential = pelectrostatic(graphite, xx, yy)
+
+		self.assertSequenceEqual(xx.shape, potential.shape)
+
+	def test_side_effects(self):
+		""" Test that mesh arrays are not written to in pelectrostatic """
+		xx, yy = np.meshgrid(np.linspace(-10, 10, 32), np.linspace(-10, 10, 32))
+		xx.setflags(write = False)
+		yy.setflags(write = False)
+		potential = pelectrostatic(graphite, xx, yy)
+
+	def test_trivial(self):
+		""" Test that the projected electrostatic potential from an empty slice of crystal is zero"""
+		xx, yy = np.meshgrid(np.linspace(-10, 10, 32), np.linspace(-10, 10, 32))
+		potential = pelectrostatic(graphite, xx, yy, bounds = (1,2))
+		self.assertTrue(np.allclose(potential, 0))
+
+if __name__ == '__main__':
+	unittest.main()