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:
u=UnitCell(a,b,c,alpha,beta,gamma)
m=numpy.linalg.inv(u.getB())
ijk=[0,0.5,0.25] # an atom (or muon!) position in the cell
xyz=numpy.dot(ijk,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 v.xyz() 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 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:
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 [np.dot(eq, 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
NN21=2*NN+1
indx=numpy.indices([NN21,NN21,NN21]).reshape((3,-1)).transpose()-NN # unit cell offsets in 3D, reduced coordinates
uc=UnitCell(5,5,6,90,90,120)
try:
sg=SpaceGroupFactory.createSpaceGroup('P 63/m')
except:
print SpaceGroupFactory.getAllSpaceGroupSymbols()
raise
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]
bigcontentsA=numpy.dot(bigcontents, uc.getBinv()) # coords in A[site, coordinate xyz]
loc0=numpy.dot([0.2,0.2,0.1], uc.getBinv()) # muon or other site of interest
d=loc0-bigcontentsA # displacement vectors
r=numpy.linalg.norm(d,axis=1)
print numpy.sort(r) # list of nearest neighbour distances