Real space crystallography

Mantid has a lot of helpful looking functions to work with space and point groups and unit cells. However all the documentation seems to be biased towards reciprocal space, which is of course useful in diffraction. With muons for example we’re much more interested in real space - such as what are the locations on Cartesian axes of atoms with reduced coordinates (x,y,z) in this and adjacent cells or what are the equivalent directions in real space of a magnetic moment oriented in some arbitrary direction in a fairly high symmetry structure.

It may be that the same maths works just the same for real and reciprocal spaces, or there could be some space groups for which there are subtle catches to beware of. So the documentation ought to be updated, or perhaps just put links to existing documentation elsewhere on the web if it uses the same notation as Mantid. And are there missing algorithms or functions?

In terms of things that are exposed to Python, there are ways to use unit cells, space- and point groups directly. With space group objects you can find equivalent positions of an atom in reduced coordinates with a script like this:

from mantid.geometry import SpaceGroupFactory

space_group = SpaceGroupFactory.createSpaceGroup('F m -3 m')
equivalents = space_group.getEquivalentPositions([0.15, 0.25, 0.11])

print equivalents

To get the reduced coordinated for the positions in cells in other directions you could add mantid.kernel.V3D-instances to each element of the list and then they could also be transformed to cartesian coordinates in Angstrom using UnitCell. The group-related things are currently documented as concept pages:

I agree that there are some pitfalls when it comes to transforming real space vs. reciprocal space coordinates and vectors, especially for everything that is hexagonal (or trigonal in hexagonal setting). Maybe we could have a discussion about what would be useful additions for your use cases.

Those algorithms work in “real detector space” whereas I’d like to work in “real crystal space” given a space group, lattice parameters and cell contents already determined. Perhaps putting “crystallography” in the title was misleading.

So far, from experimentation. I can do:
ijk=[0,0.5,0.25] # an atom (or muon!) position in the cell,m) # coords in Angstroms

The documentation in “Concepts - Lattice” hints at the existence of a “u.getBinv()” method but that doesn’t seem to exist in reality, hence the explicit inversion.

p.s. Mantid’s “V3D” type is fiddly to convert to a plain 3 element list or array. Should there be a method that returns a three element list, or even that it behaves as such if used in a list context such as nv=numpy.array(v) ?

pps. Is there an algorithm or method to verify that a (user supplied) unit cell and space group are compatible (e.g. angles 90 or 120 degrees, a=b, etc) before my routine goes ahead and produces junk output?

While trying out the script I posted before I also noticed that V3D is a bit hard to convert. On the C++ level it already has a method to convert to a std::vector, I think it should be possible to expose that to Python so that a numpy array is returned.

I opened an issue for this: Export std::vector conversion from Kernel::V3D to Python · Issue #16160 · mantidproject/mantid · GitHub

I looked into the V3D issue this afternoon and it turned out that it was possible to make it behave as you suggested.

As soon as the pull request is checked and merged the behaviour you requested will be available in the nightly build and then also in the upcoming release:

The above pull request has made it into the master branch, the new V3D-bahvior will be part of the upcoming release.

With respect to checking space group and unit cell compatibility, there is now a second pull request where this functionality has been implemented.

After these changes you can simply do this:

cell = UnitCell(5, 5, 6, 90, 90, 120)

sg_cubic = SpaceGroupFactory.createSpaceGroup("P m -3 m")
sg_hexagonal = SpaceGroupFactory.createSpaceGroup("P 6/m")

print "Metric compatible with cubic symmetry:", sg_cubic.isAllowedUnitCell(cell)
print "Metric compatible with hexagonal symmetry:", sg_hexagonal.isAllowedUnitCell(cell)

The documentation in the point and space group concept page is also updated as part of that pull request to describe this feature.

I hope this matches what you had in mind.

To make this small series of fixes complete, adding the Python export for UnitCell::getBinv:

With that last change it is possible to write a function like this:

import numpy as np

def get_atom_coordinates_in_angstrom(position, space_group, cell):
    if not space_group.isAllowedUnitCell(cell):
        raise RuntimeError('The supplied unit cell is not compatible with space group {0}.'.format(space_group.getHMSymbol()))
    return [, cell.getBinv()) for eq in space_group.getEquivalentPositions(position)]
sg = SpaceGroupFactory.createSpaceGroup('P m m m')
cell = UnitCell(4, 6, 8, 90, 90, 90)

print get_atom_coordinates_in_angstrom([0.0, 0.25, 0.125], sg, cell)

It takes a position in reduced coordinates and returns a list of all equivalent coordinates in the cell in Angstrom, using the procedure you outlined above.

Given that release 3.7 is very close now, we could add more functionality (better support for working with real space directions, etc.) in the next release 3.8. If you have any suggestions I’d be happy to hear about them.

Thanks. That does what I need. It’s now easy to process a lot of sites in one go with a NumPy array:

import numpy

NN=1 # nearest neighbour cells to consider in each direction
indx=numpy.indices([NN21,NN21,NN21]).reshape((3,-1)).transpose()-NN # unit cell offsets in 3D, reduced coordinates

	sg=SpaceGroupFactory.createSpaceGroup('P 63/m')
	print SpaceGroupFactory.getAllSpaceGroupSymbols()
if not sg.isAllowedUnitCell(uc):
	raise RuntimeError('The supplied unit cell is not compatible with space group {0}.'.format(sg.getHMSymbol()))

contents2=numpy.array(sg.getEquivalentPositions([0.25,0,0.25])) # one unit cell

bigcontents=(contents2[:,numpy.newaxis,:]+indx[numpy.newaxis,:,:]).reshape((-1,3)) # many unit cells, reduced coords[site, coordinate xyz], uc.getBinv()) # coords in A[site, coordinate xyz][0.2,0.2,0.1], uc.getBinv()) # muon or other site of interest

d=loc0-bigcontentsA # displacement vectors
print numpy.sort(r) # list of nearest neighbour distances