Commit 565e4c42 authored by Lynch, Vickie's avatar Lynch, Vickie
Browse files

Merge pull request #14971 from mantidproject/11268_create_crystal_structure_loader

Create algorithm to load crystal structure from CIF
parents fc97a8c5 5c028c05
......@@ -14,23 +14,45 @@
namespace Mantid {
namespace Geometry {
/// Functor for fuzzy comparison of V3D-objects using Kernel::Tolerance
struct MANTID_GEOMETRY_DLL FuzzyV3DLessThan {
bool operator()(const Kernel::V3D &lhs, const Kernel::V3D &rhs) {
if (fabs(lhs.X() - rhs.X()) > Kernel::Tolerance) {
/// Equality-functor for comparison of atom positions with specifiable precision
class MANTID_GEOMETRY_DLL AtomPositionsEqual {
public:
AtomPositionsEqual(double precision = 1.e-4) : m_precision(precision) {}
bool operator()(const Kernel::V3D &lhs, const Kernel::V3D &rhs) const {
return !(fabs(lhs.X() - rhs.X()) > m_precision ||
fabs(lhs.Y() - rhs.Y()) > m_precision ||
fabs(lhs.Z() - rhs.Z()) > m_precision);
}
private:
double m_precision;
};
/// Less-than-functor for comparison of atom positions with specifiable
/// precision
class MANTID_GEOMETRY_DLL AtomPositionsLessThan {
public:
AtomPositionsLessThan(double precision = 1.e-4) : m_precision(precision) {}
bool operator()(const Kernel::V3D &lhs, const Kernel::V3D &rhs) const {
if (fabs(lhs.X() - rhs.X()) > m_precision) {
return lhs.X() < rhs.X();
}
if (fabs(lhs.Y() - rhs.Y()) > Kernel::Tolerance) {
if (fabs(lhs.Y() - rhs.Y()) > m_precision) {
return lhs.Y() < rhs.Y();
}
if (fabs(lhs.Z() - rhs.Z()) > Kernel::Tolerance) {
if (fabs(lhs.Z() - rhs.Z()) > m_precision) {
return lhs.Z() < rhs.Z();
}
return false;
}
private:
double m_precision;
};
/**
......
......@@ -81,8 +81,9 @@ public:
}
// Use fuzzy compare with the same condition as V3D::operator==().
std::sort(equivalents.begin(), equivalents.end(), FuzzyV3DLessThan());
equivalents.erase(std::unique(equivalents.begin(), equivalents.end()),
std::sort(equivalents.begin(), equivalents.end(), AtomPositionsLessThan());
equivalents.erase(std::unique(equivalents.begin(), equivalents.end(),
AtomPositionsEqual()),
equivalents.end());
return equivalents;
......
......@@ -91,7 +91,7 @@ std::vector<Kernel::V3D> Group::operator*(const Kernel::V3D &vector) const {
result.push_back(Geometry::getWrappedVector((*op) * vector));
}
std::sort(result.begin(), result.end(), FuzzyV3DLessThan());
std::sort(result.begin(), result.end(), AtomPositionsLessThan());
result.erase(std::unique(result.begin(), result.end()), result.end());
return result;
......
......@@ -105,9 +105,11 @@ Group_const_sptr SpaceGroup::getSiteSymmetryGroup(const V3D &position) const {
std::vector<SymmetryOperation> siteSymmetryOps;
AtomPositionsEqual comparator;
for (auto op = m_allOperations.begin(); op != m_allOperations.end(); ++op) {
if (Geometry::getWrappedVector((*op) * wrappedPosition) ==
wrappedPosition) {
if (comparator(Geometry::getWrappedVector((*op) * wrappedPosition),
wrappedPosition)) {
siteSymmetryOps.push_back(*op);
}
}
......
......@@ -185,8 +185,8 @@ public:
TS_ASSERT_EQUALS(two.getCoordinateSystem(), Group::Hexagonal);
}
void testFuzzyV3DLessThan() {
FuzzyV3DLessThan lessThan;
void testAtomPositionLessThan() {
AtomPositionsLessThan lessThan(1.e-6);
V3D v1(0.654321, 0.0, 0.0);
V3D v2(0.654320, 0.0, 0.0);
......@@ -220,6 +220,33 @@ public:
TS_ASSERT(!lessThan(v7, v1));
}
void testAtomPositionEqual() {
AtomPositionsEqual equal(1.e-6);
V3D v1(0.654321, 0.0, 0.0);
V3D v2(0.654320, 0.0, 0.0);
TS_ASSERT(!equal(v2, v1));
// 7th digit is not compared.
V3D v3(0.6543211, 0.0, 0.0);
TS_ASSERT(equal(v1, v3));
TS_ASSERT(equal(v3, v1));
// Same for y
V3D v4(0.654321, 0.0000010001, 0.0);
TS_ASSERT(!equal(v1, v4));
V3D v5(0.654321, 0.0000001, 0.0);
TS_ASSERT(equal(v5, v1));
// Same for z
V3D v6(0.654321, 0.0, 0.0000010001);
TS_ASSERT(!equal(v1, v6));
V3D v7(0.654321, 0.0, 0.0000001);
TS_ASSERT(equal(v1, v7));
}
void testContainsOperation() {
Group group("x,y,z; -x,-y,-z");
......
# pylint: disable=no-init,too-few-public-methods
from mantid.kernel import *
from mantid.simpleapi import *
from mantid.api import *
from mantid.geometry import SpaceGroupFactory, CrystalStructure
import re
# pylint: disable=invalid-name
def removeErrorEstimateFromNumber(numberString):
errorBegin = numberString.find('(')
if errorBegin == -1:
return numberString
return numberString[:errorBegin]
class CrystalStructureBuilder(object):
'''
Helper class that builds a CrystalStructure file from the result
of the ReadCif-function provided by PyCifRW.
For testing purposes, dictionaries with the appropriate data can
be passed in as well, so the source of the parsed data is replaceable.
'''
def __init__(self, cifFile=None):
if cifFile is not None:
cifData = cifFile[cifFile.keys()[0]]
self.spaceGroup = self._getSpaceGroup(cifData)
self.unitCell = self._getUnitCell(cifData)
self.atoms = self._getAtoms(cifData)
def getCrystalStructure(self):
return CrystalStructure(self.unitCell, self.spaceGroup, self.atoms)
def _getSpaceGroup(self, cifData):
try:
return self._getSpaceGroupFromString(cifData)
# pylint: disable=unused-variable
except (RuntimeError, ValueError) as error:
try:
return self._getSpaceGroupFromNumber(cifData)
# pylint: disable=unused-variable,invalid-name
except RuntimeError as e:
raise RuntimeError(
'Can not create space group from supplied CIF-file. You could try to modify the HM-symbol'\
'to contain spaces between the components.\n'\
'Keys to look for: _space_group_name_H-M_alt, _symmetry_space_group_name_H-M')
def _getSpaceGroupFromString(self, cifData):
# Try two possibilities for space group symbol. If neither is present, throw a RuntimeError.
rawSpaceGroupSymbol = [str(cifData[x]) for x in
[u'_space_group_name_h-m_alt', u'_symmetry_space_group_name_h-m'] if
x in cifData.keys()]
if len(rawSpaceGroupSymbol) == 0:
raise RuntimeError('No space group symbol in CIF.')
cleanSpaceGroupSymbol = self._getCleanSpaceGroupSymbol(rawSpaceGroupSymbol[0])
# If the symbol is not registered, throw as well.
return SpaceGroupFactory.createSpaceGroup(cleanSpaceGroupSymbol).getHMSymbol()
def _getCleanSpaceGroupSymbol(self, rawSpaceGroupSymbol):
# Remove :1 and :H from the symbol. Those are not required at the moment because they are the default.
removalRe = re.compile(':[1H]', re.IGNORECASE)
return re.sub(removalRe, '', rawSpaceGroupSymbol).strip()
def _getSpaceGroupFromNumber(self, cifData):
spaceGroupNumber = [int(cifData[x]) for x in
[u'_space_group_it_number', u'_symmetry_int_tables_number'] if
x in cifData.keys()]
if len(spaceGroupNumber) == 0:
raise RuntimeError('No space group symbol in CIF.')
possibleSpaceGroupSymbols = SpaceGroupFactory.subscribedSpaceGroupSymbols(spaceGroupNumber[0])
if len(possibleSpaceGroupSymbols) != 1:
raise RuntimeError(
'Can not use space group number to determine space group for no. {0}'.format(spaceGroupNumber))
return SpaceGroupFactory.createSpaceGroup(possibleSpaceGroupSymbols[0]).getHMSymbol()
def _getUnitCell(self, cifData):
unitCellComponents = [u'_cell_length_a', u'_cell_length_b', u'_cell_length_c',
u'_cell_angle_alpha', u'_cell_angle_beta', u'_cell_angle_gamma']
unitCellValueMap = dict([(str(x), removeErrorEstimateFromNumber(str(cifData[x]))) if x in cifData.keys()
else (str(x), None) for x in
unitCellComponents])
if unitCellValueMap['_cell_length_a'] is None:
raise RuntimeError('The a-parameter of the unit cell is not specified in the supplied CIF.\n' \
'Key to look for: _cell_length_a')
replacementMap = {
'_cell_length_b': str(unitCellValueMap['_cell_length_a']),
'_cell_length_c': str(unitCellValueMap['_cell_length_a']),
'_cell_angle_alpha': '90.0', '_cell_angle_beta': '90.0', '_cell_angle_gamma': '90.0'}
unitCellValues = [
unitCellValueMap[str(key)] if unitCellValueMap[str(key)] is not None else replacementMap[str(key)] for key
in unitCellComponents]
return ' '.join(unitCellValues)
def _getAtoms(self, cifData):
atomFieldsRequirements = [(u'_atom_site_label', True),
(u'_atom_site_fract_x', True),
(u'_atom_site_fract_y', True),
(u'_atom_site_fract_z', True),
(u'_atom_site_occupancy', False),
(u'_atom_site_U_iso_or_equiv', False)]
atomFields = []
for field, required in atomFieldsRequirements:
if field in cifData.keys():
atomFields.append(field)
else:
if required:
raise RuntimeError(
'Mandatory field {0} not found in CIF-file.' \
'Please check the atomic position definitions.'.format(field))
atomLists = [cifData[x] for x in atomFields]
atomLines = []
for atomLine in zip(*atomLists):
stringAtomLine = [str(x) for x in atomLine]
cleanLine = [self._getCleanAtomSymbol(stringAtomLine[0])] + [removeErrorEstimateFromNumber(x) for x in
list(stringAtomLine[1:])]
atomLines.append(' '.join(cleanLine))
return ';'.join(atomLines)
def _getCleanAtomSymbol(self, atomSymbol):
nonCharacterRe = re.compile('[^a-z]', re.IGNORECASE)
return re.sub(nonCharacterRe, '', atomSymbol)
class UBMatrixBuilder(object):
def __init__(self, cifFile=None):
if cifFile is not None:
cifData = cifFile[cifFile.keys()[0]]
self._ubMatrix = self._getUBMatrix(cifData)
def getUBMatrix(self):
return self._ubMatrix
def _getUBMatrix(self, cifData):
ubMatrixKeys = [u'_diffrn_orient_matrix_ub_11', u'_diffrn_orient_matrix_ub_12', u'_diffrn_orient_matrix_ub_13',
u'_diffrn_orient_matrix_ub_21', u'_diffrn_orient_matrix_ub_22', u'_diffrn_orient_matrix_ub_23',
u'_diffrn_orient_matrix_ub_31', u'_diffrn_orient_matrix_ub_32', u'_diffrn_orient_matrix_ub_33']
ubValues = [str(cifData[key]) if key in cifData.keys() else None for key in ubMatrixKeys]
if None in ubValues:
raise RuntimeError('Can not load UB matrix from CIF, values are missing.')
return ','.join(ubValues)
class LoadCIF(PythonAlgorithm):
def category(self):
return "Diffraction\\DataHandling"
def name(self):
return "LoadCIF"
def summary(self):
return "This algorithm loads a CIF file using the PyCifRW package and assigns a CrystalStructure to the sample of the workspace."
def PyInit(self):
self.declareProperty(
WorkspaceProperty(name='Workspace',
defaultValue='', direction=Direction.InOut),
doc='Workspace into which the crystal structure is placed.')
self.declareProperty(
FileProperty(name='InputFile',
defaultValue='',
action=FileAction.Load,
extensions=['cif']),
doc='A CIF file containing a crystal structure.')
self.declareProperty('LoadUBMatrix', False,
direction=Direction.Input,
doc='Load UB-matrix from CIF file if available.')
def PyExec(self):
cifFileName = self.getProperty('InputFile').value
workspace = self.getProperty('Workspace').value
# Try to parse cif file using PyCifRW
parsedCifFile = ReadCif(cifFileName)
self._setCrystalStructureFromCifFile(workspace, parsedCifFile)
ubOption = self.getProperty('LoadUBMatrix').value
if ubOption:
self._setUBMatrixFromCifFile(workspace, parsedCifFile)
def _setCrystalStructureFromCifFile(self, workspace, cifFile):
crystalStructure = self._getCrystalStructureFromCifFile(cifFile)
workspace.sample().setCrystalStructure(crystalStructure)
def _getCrystalStructureFromCifFile(self, cifFile):
builder = CrystalStructureBuilder(cifFile)
crystalStructure = builder.getCrystalStructure()
self.log().information('''Loaded the following crystal structure:
Unit cell: {0}
Space group: {1}
Atoms: {2}
'''.format(builder.unitCell, builder.spaceGroup, builder.atoms))
return crystalStructure
def _setUBMatrixFromCifFile(self, workspace, cifFile):
ubMatrix = self._getUBMatrixFromCifFile(cifFile)
setUBAlgorithm = self.createChildAlgorithm('SetUB')
setUBAlgorithm.setProperty('Workspace', workspace)
setUBAlgorithm.setProperty('UB', ubMatrix)
setUBAlgorithm.execute()
def _getUBMatrixFromCifFile(self, cifFile):
builder = UBMatrixBuilder(cifFile)
return builder.getUBMatrix()
try:
from CifFile import ReadCif
AlgorithmFactory.subscribe(LoadCIF)
except ImportError:
logger.debug('Failed to subscribe algorithm LoadCIF; Python package PyCifRW' \
'may be missing (https://pypi.python.org/pypi/PyCifRW/4.1)')
......@@ -92,6 +92,7 @@ set ( TEST_PY_FILES
PoldiMergeTest.py
VesuvioResolutionTest.py
PoldiCreatePeaksFromFileTest.py
LoadCIFTest.py
)
check_tests_valid ( ${CMAKE_CURRENT_SOURCE_DIR} ${TEST_PY_FILES} )
......
# pylint: disable=no-init,too-many-public-methods,invalid-name,protected-access
import unittest
from testhelpers import assertRaisesNothing
from LoadCIF import UBMatrixBuilder, CrystalStructureBuilder
from mantid.api import AlgorithmFactory
def merge_dicts(lhs, rhs):
merged = lhs.copy()
merged.update(rhs)
return merged
class CrystalStructureBuilderTestSpaceGroup(unittest.TestCase):
def setUp(self):
self.builder = CrystalStructureBuilder()
def test_getSpaceGroupFromString_valid_no_exceptions(self):
valid_new = {u'_space_group_name_h-m_alt': u'P m -3 m'}
valid_old = {u'_symmetry_space_group_name_h-m': u'P m -3 m'}
assertRaisesNothing(self, self.builder._getSpaceGroupFromString, cifData=valid_old)
assertRaisesNothing(self, self.builder._getSpaceGroupFromString, cifData=valid_new)
assertRaisesNothing(self, self.builder._getSpaceGroupFromString, cifData=merge_dicts(valid_new, valid_old))
def test_getSpaceGroupFromString_valid_correct_value(self):
valid_new = {u'_space_group_name_h-m_alt': u'P m -3 m'}
valid_old = {u'_symmetry_space_group_name_h-m': u'P m -3 m'}
valid_old_different = {u'_symmetry_space_group_name_h-m': u'F d d d'}
invalid_old = {u'_symmetry_space_group_name_h-m': u'invalid'}
self.assertEqual(self.builder._getSpaceGroupFromString(valid_new), 'P m -3 m')
self.assertEqual(self.builder._getSpaceGroupFromString(valid_old), 'P m -3 m')
self.assertEqual(self.builder._getSpaceGroupFromString(merge_dicts(valid_new, valid_old)), 'P m -3 m')
self.assertEqual(self.builder._getSpaceGroupFromString(merge_dicts(valid_new, valid_old_different)), 'P m -3 m')
self.assertEqual(self.builder._getSpaceGroupFromString(merge_dicts(valid_new, invalid_old)), 'P m -3 m')
def test_getSpaceGroupFromString_invalid(self):
valid_old = {u'_symmetry_space_group_name_h-m': u'P m -3 m'}
invalid_new = {u'_space_group_name_h-m_alt': u'invalid'}
invalid_old = {u'_symmetry_space_group_name_h-m': u'invalid'}
self.assertRaises(RuntimeError, self.builder._getSpaceGroupFromString, cifData={})
self.assertRaises(ValueError, self.builder._getSpaceGroupFromString, cifData=invalid_new)
self.assertRaises(ValueError, self.builder._getSpaceGroupFromString, cifData=invalid_old)
self.assertRaises(ValueError, self.builder._getSpaceGroupFromString,
cifData=merge_dicts(invalid_new, valid_old))
def test_getCleanSpaceGroupSymbol(self):
fn = self.builder._getCleanSpaceGroupSymbol
self.assertEqual(fn('P m -3 m :1'), 'P m -3 m')
self.assertEqual(fn('P m -3 m :H'), 'P m -3 m')
def test_getSpaceGroupFromNumber_invalid(self):
invalid_old = {u'_symmetry_int_tables_number': u'400'}
invalid_new = {u'_space_group_it_number': u'400'}
self.assertRaises(RuntimeError, self.builder._getSpaceGroupFromNumber, cifData={})
self.assertRaises(RuntimeError, self.builder._getSpaceGroupFromNumber, cifData=invalid_old)
self.assertRaises(RuntimeError, self.builder._getSpaceGroupFromNumber, cifData=invalid_new)
class CrystalStructureBuilderTestUnitCell(unittest.TestCase):
def setUp(self):
self.builder = CrystalStructureBuilder()
def test_getUnitCell_invalid(self):
invalid_no_a = {u'_cell_length_b': u'5.6'}
self.assertRaises(RuntimeError, self.builder._getUnitCell, cifData={})
self.assertRaises(RuntimeError, self.builder._getUnitCell, cifData=invalid_no_a)
def test_getUnitCell_cubic(self):
cell = {u'_cell_length_a': u'5.6'}
self.assertEqual(self.builder._getUnitCell(cell), '5.6 5.6 5.6 90.0 90.0 90.0')
def test_getUnitCell_tetragonal(self):
cell = {u'_cell_length_a': u'5.6', u'_cell_length_c': u'2.3'}
self.assertEqual(self.builder._getUnitCell(cell), '5.6 5.6 2.3 90.0 90.0 90.0')
def test_getUnitCell_orthorhombic(self):
cell = {u'_cell_length_a': u'5.6', u'_cell_length_b': u'1.6', u'_cell_length_c': u'2.3'}
self.assertEqual(self.builder._getUnitCell(cell), '5.6 1.6 2.3 90.0 90.0 90.0')
def test_getUnitCell_hexagonal(self):
cell = {u'_cell_length_a': u'5.6', u'_cell_length_c': u'2.3', u'_cell_angle_gamma': u'120.0'}
cell_errors = {u'_cell_length_a': u'5.6(1)', u'_cell_length_c': u'2.3(1)', u'_cell_angle_gamma': u'120.0'}
self.assertEqual(self.builder._getUnitCell(cell), '5.6 5.6 2.3 90.0 90.0 120.0')
self.assertEqual(self.builder._getUnitCell(cell_errors), '5.6 5.6 2.3 90.0 90.0 120.0')
class CrystalStructureBuilderTestAtoms(unittest.TestCase):
def setUp(self):
self.builder = CrystalStructureBuilder()
def test_getAtoms_required_keys(self):
mandatoryKeys = dict([(u'_atom_site_label', [u'Si']),
(u'_atom_site_fract_x', [u'1/8']),
(u'_atom_site_fract_y', [u'1/8']),
(u'_atom_site_fract_z', [u'1/8'])])
for key in mandatoryKeys:
tmp = mandatoryKeys.copy()
del tmp[key]
self.assertRaises(RuntimeError, self.builder._getAtoms, cifData=tmp)
def test_getAtoms_correct(self):
data = dict([(u'_atom_site_label', [u'Si', u'Al']),
(u'_atom_site_fract_x', [u'1/8', u'0.34(1)']),
(u'_atom_site_fract_y', [u'1/8', u'0.56(2)']),
(u'_atom_site_fract_z', [u'1/8', u'0.23(2)']),
(u'_atom_site_occupancy', [u'1.0', u'1.0(0)']),
(u'_atom_site_U_iso_or_equiv', [u'0.01', u'0.02'])])
self.assertEqual(self.builder._getAtoms(data), 'Si 1/8 1/8 1/8 1.0 0.01;Al 0.34 0.56 0.23 1.0 0.02')
class UBMatrixBuilderTest(unittest.TestCase):
def setUp(self):
self.builder = UBMatrixBuilder()
self.valid_matrix = {u'_diffrn_orient_matrix_ub_11': u'-0.03',
u'_diffrn_orient_matrix_ub_12': u'0.13',
u'_diffrn_orient_matrix_ub_13': u'0.31',
u'_diffrn_orient_matrix_ub_21': u'0.01',
u'_diffrn_orient_matrix_ub_22': u'-0.31',
u'_diffrn_orient_matrix_ub_23': u'0.14',
u'_diffrn_orient_matrix_ub_31': u'0.34',
u'_diffrn_orient_matrix_ub_32': u'0.02',
u'_diffrn_orient_matrix_ub_33': u'0.02'}
def test_getUBMatrix_invalid(self):
for key in self.valid_matrix:
tmp = self.valid_matrix.copy()
del tmp[key]
self.assertRaises(RuntimeError, self.builder._getUBMatrix, cifData=tmp)
def test_getUBMatrix_correct(self):
self.assertEqual(self.builder._getUBMatrix(self.valid_matrix), '-0.03,0.13,0.31,0.01,-0.31,0.14,0.34,0.02,0.02')
if __name__ == '__main__':
# Only test if algorithm is registered (PyCifRW dependency).
if AlgorithmFactory.exists("LoadCIF"):
unittest.main()
.. algorithm::
.. summary::
.. alias::
.. properties::
Description
-----------
This algorithm uses ``PyCifRW`` to parse the supplied CIF-file and extract the information necessary to construct a
``CrystalStructure`` object that is attached to the supplied workspace. For successfully loading a crystal structure,
the CIF-file must contain some mandatory fields, describing the three components that define the structure.
**Unit cell**
For the unit cell, at least *_cell_length_a* must be present. If *_cell_length_b* or *_cell_length_c* are absent, the
value is replaced with the one of *_cell_length_a*. Any absent item of *_cell_angle_alpha*, *_cell_angle_beta* or
*_cell_angle_gamma* is missing, it's replaced with 90.
**Space group**
The space group is currently loaded based on one of the fields *_space_group_name_H-M_alt* or
*_symmetry_space_group_name_H-M*, where the former has precedence over the latter. The value has to correspond to one
of the registered space groups in Mantid (see the corresponding :ref:`concept page <Point and space groups>`).
If neither of the fields is present, the loader tries to recover by using the space group number in either
*_symmetry_int_tables_number* or *_space_group_it_number* where again, the former has precedence. This only works if
there is only one setting registered for the specified space group type, because otherwise the algorithm can not decide
which one to use.