diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/BasicHKLFilters.h b/Framework/Geometry/inc/MantidGeometry/Crystal/BasicHKLFilters.h new file mode 100644 index 0000000000000000000000000000000000000000..8081a3ebe33d623b95637c87a4995dfab5e02f74 --- /dev/null +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/BasicHKLFilters.h @@ -0,0 +1,162 @@ +#ifndef MANTID_GEOMETRY_BASICHKLFILTERS_H_ +#define MANTID_GEOMETRY_BASICHKLFILTERS_H_ + +#include "MantidGeometry/DllConfig.h" +#include "MantidGeometry/Crystal/HKLFilter.h" +#include "MantidGeometry/Crystal/ReflectionCondition.h" +#include "MantidGeometry/Crystal/StructureFactorCalculator.h" +#include "MantidGeometry/Crystal/SpaceGroup.h" +#include "MantidGeometry/Crystal/UnitCell.h" + +#include <strstream> + +namespace Mantid { +namespace Geometry { + +/** BasicHKLFilters + + This file contains some implementations of HKLFilter that + provide filtering based on things like d-value, space group + or centering. + + A common use would be to generate a specific list of HKLs, + for example all reflections that are allowed according to a certain + range of d-values and the reflection conditions of a space group: + + HKLFilter_const_sptr filter = + boost::make_shared<HKLFilterDRange>(unitCell, 0.5) + & boost::make_shared<HKLFilterSpaceGroup>(spaceGroup); + + HKLGenerator gen(unitCell, 0.5); + std::vector<V3D> hkls; + + std::remove_copy_if(gen.begin(), gen.end(), std::back_inserter(hkls), + (~filter)->fn()); + + An existing list of HKLs could be checked for indices that match the + reflection conditions of a space group: + + HKLFilter_const_sptr sgFilter = + boost::make_shared<HKLFilterSpaceGroup>(spaceGroup); + + auto matchingHKLCount = std::count_if(hkls.begin(), hkls.end(), + sgFilter->fn()); + + auto violatingHKLCount = std::count_if(hkls.begin(), hkls.end(), + (~sgFilter)->fn()); + + Combining HKLGenerator and different HKLFilters provides a very flexible + system for creating and processing specific sets of Miller indices that + is easy to expand by adding other HKLFilters. + + @author Michael Wedel, ESS + @date 23/09/2015 + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> +*/ + +/** + * A class to filter HKLs by their d-values + * + * This class takes a UnitCell object and calculates the spacing of + * the lattice planes for each HKL. If the lattice spacing is within + * the spcified range of values, the reflection is allowed. + * + * If the first constructor with only dMin is used, dMax is taken to + * be the lattice parameter with the largest value. There can not be + * a greater interplanar spacing than that value. + */ +class MANTID_GEOMETRY_DLL HKLFilterDRange : public HKLFilter { +public: + HKLFilterDRange(const UnitCell &cell, double dMin); + HKLFilterDRange(const UnitCell &cell, double dMin, double dMax); + + std::string getDescription() const; + bool isAllowed(const Kernel::V3D &hkl) const; + +private: + void checkProperDRangeValues(); + + UnitCell m_cell; + double m_dmin, m_dmax; +}; + +/** + * A class to filter HKLs according to a space group + * + * HKLFilterSpaceGroup stores a SpaceGroup object and marks those + * reflections as allowed that are allowed according to the + * reflection conditions of the space group. + */ +class MANTID_GEOMETRY_DLL HKLFilterSpaceGroup : public HKLFilter { +public: + HKLFilterSpaceGroup(const SpaceGroup_const_sptr &spaceGroup); + + std::string getDescription() const; + bool isAllowed(const Kernel::V3D &hkl) const; + +protected: + SpaceGroup_const_sptr m_spaceGroup; +}; + +/** + * A class to filter HKLs according to structure factor magnitudes + * + * This filter takes a StructureFactorCalculator and calculates the + * structure factor for each HKL. If F^2 is larger than the specified + * minimum, the reflection is considered allowed. The default minimum + * is 1e-6. + */ +class MANTID_GEOMETRY_DLL HKLFilterStructureFactor : public HKLFilter { +public: + HKLFilterStructureFactor(const StructureFactorCalculator_sptr &calculator, + double fSquaredMin = 1.0e-6); + + std::string getDescription() const; + bool isAllowed(const Kernel::V3D &hkl) const; + +protected: + StructureFactorCalculator_sptr m_calculator; + double m_fSquaredMin; +}; + +/** + * A class to filter HKLs according to a lattice centering + * + * HKLFilterCentering is a filter that stores a ReflectionCondition object + * internally and filters the HKLs according to that. + */ +class MANTID_GEOMETRY_DLL HKLFilterCentering : public HKLFilter { +public: + HKLFilterCentering(const ReflectionCondition_sptr ¢ering); + + std::string getDescription() const; + bool isAllowed(const Kernel::V3D &hkl) const; + +protected: + ReflectionCondition_sptr m_centering; +}; + +} // namespace Geometry +} // namespace Mantid + +#endif /* MANTID_GEOMETRY_BASICHKLFILTERS_H_ */ diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/HKLFilter.h b/Framework/Geometry/inc/MantidGeometry/Crystal/HKLFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..ca812aebdefeb7fdaa31219407f718d6becdfa2b --- /dev/null +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/HKLFilter.h @@ -0,0 +1,175 @@ +#ifndef MANTID_GEOMETRY_HKLFILTER_H_ +#define MANTID_GEOMETRY_HKLFILTER_H_ + +#include "MantidGeometry/DllConfig.h" +#include "MantidKernel/V3D.h" + +#include <boost/shared_ptr.hpp> +#include <functional> + +namespace Mantid { +namespace Geometry { + +/** HKLFilter + + There are many ways to filter lists of Miller indices HKL. In order + to be able to use HKLGenerator with arbitrary filters, HKLFilter provides + a general interface for such filters. + + The abstract base class HKLFilter defines a pure virtual method + HKLFilter::isAllowed(), which takes a V3D as argument and returns a boolean. + Implementing classes can then implement this method, wrapping very different + concepts of checking HKLs for certain characteristics. + + There are two general ways of using HKLFilters. When used "standalone", + the isAllowed()-function can be used directly: + + if(filter->isAllowed(hkl)) { + // do something + } + + For interoperability with STL-algorithms, HKLFilter provides a method that + returns a function-object with the filter: + + std::copy_if(generator.begin(), generator.end(), + hkls.begin(), filter->fn()); + + Often, it's not enough to filter a list by one criterion. To this end, + there are two implementations of HKLFilter that provide binary logic + operations + "and" (HKLFilterAnd) and "or" (HKLFilterOr). They can be constructed from two + HKLFilter_const_sptrs, or, more conveniently, by using the operators & and | + directly with those types. This makes it possible to combine filters in many + ways: + + HKLFilter_const_sptr filter = (filter1 | filter2) & filter3; + + Lastly, the unary logic operation "not" (HKLFilterNot) is implemented. This is + important for usability with the pre-C++11 algorithm std::remove_copy_if that + works logically reversed compared to std::copy_if: + + std::remove_copy_if(generator.begin(), generator.end(), hkls.begin(), + (~filter)->fn()); + + For actual implementations of HKLFilter, check out BasicHKLFilters, where some + important filters are defined. + + @author Michael Wedel, ESS + @date 06/09/2015 + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> +*/ +class MANTID_GEOMETRY_DLL HKLFilter { +public: + virtual ~HKLFilter() {} + + std::function<bool(const Kernel::V3D &)> fn() const; + + virtual std::string getDescription() const = 0; + virtual bool isAllowed(const Kernel::V3D &hkl) const = 0; +}; + +typedef boost::shared_ptr<const HKLFilter> HKLFilter_const_sptr; + +/// Base class for unary logic operations for HKLFilter. +class MANTID_GEOMETRY_DLL HKLFilterUnaryLogicOperation : public HKLFilter { +public: + HKLFilterUnaryLogicOperation(const HKLFilter_const_sptr &filter); + ~HKLFilterUnaryLogicOperation() {} + + /// Returns the operand of the function. + const HKLFilter_const_sptr &getOperand() const { return m_operand; } + +protected: + HKLFilter_const_sptr m_operand; +}; + +/// Logical "Not"-operation for HKLFilter. +class MANTID_GEOMETRY_DLL HKLFilterNot : public HKLFilterUnaryLogicOperation { +public: + /// Constructor, calls base class constructor, throws exception if filter is a + /// null pointer. + HKLFilterNot(const HKLFilter_const_sptr &filter) + : HKLFilterUnaryLogicOperation(filter) {} + ~HKLFilterNot() {} + + std::string getDescription() const; + bool isAllowed(const Kernel::V3D &hkl) const; +}; + +/// Base class for binary logic operations for HKLFilter. +class MANTID_GEOMETRY_DLL HKLFilterBinaryLogicOperation : public HKLFilter { +public: + HKLFilterBinaryLogicOperation(const HKLFilter_const_sptr &lhs, + const HKLFilter_const_sptr &rhs); + virtual ~HKLFilterBinaryLogicOperation() {} + + /// Returns the left-hand side operand of the operation. + const HKLFilter_const_sptr &getLHS() const { return m_lhs; } + + /// Returns the right-hand side operand of the operation. + const HKLFilter_const_sptr &getRHS() const { return m_rhs; } + +protected: + HKLFilter_const_sptr m_lhs; + HKLFilter_const_sptr m_rhs; +}; + +/// Logical "And"-operation for HKLFilter. +class MANTID_GEOMETRY_DLL HKLFilterAnd : public HKLFilterBinaryLogicOperation { +public: + /// Constructor, calls base class constructor, throws exception if either of + /// the operands is null. + HKLFilterAnd(const HKLFilter_const_sptr &lhs, const HKLFilter_const_sptr &rhs) + : HKLFilterBinaryLogicOperation(lhs, rhs) {} + ~HKLFilterAnd() {} + + std::string getDescription() const; + bool isAllowed(const Kernel::V3D &hkl) const; +}; + +/// Logical "Or"-operation for HKLFilter. +class MANTID_GEOMETRY_DLL HKLFilterOr : public HKLFilterBinaryLogicOperation { +public: + /// Constructor, calls base class constructor, throws exception if either of + /// the operands is null. + HKLFilterOr(const HKLFilter_const_sptr &lhs, const HKLFilter_const_sptr &rhs) + : HKLFilterBinaryLogicOperation(lhs, rhs) {} + ~HKLFilterOr() {} + + std::string getDescription() const; + bool isAllowed(const Kernel::V3D &hkl) const; +}; + +MANTID_GEOMETRY_DLL const HKLFilter_const_sptr +operator~(const HKLFilter_const_sptr &filter); + +MANTID_GEOMETRY_DLL const HKLFilter_const_sptr +operator&(const HKLFilter_const_sptr &lhs, const HKLFilter_const_sptr &rhs); + +MANTID_GEOMETRY_DLL const HKLFilter_const_sptr +operator|(const HKLFilter_const_sptr &lhs, const HKLFilter_const_sptr &rhs); + +} // namespace Geometry +} // namespace Mantid + +#endif /* MANTID_GEOMETRY_HKLFILTER_H_ */ diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/HKLGenerator.h b/Framework/Geometry/inc/MantidGeometry/Crystal/HKLGenerator.h new file mode 100644 index 0000000000000000000000000000000000000000..5c4259677314690be059cfb72b1d252d8f5e19c8 --- /dev/null +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/HKLGenerator.h @@ -0,0 +1,182 @@ +#ifndef MANTID_GEOMETRY_HKLGENERATOR_H_ +#define MANTID_GEOMETRY_HKLGENERATOR_H_ + +#include "MantidGeometry/DllConfig.h" +#include "MantidGeometry/Crystal/UnitCell.h" +#include "MantidKernel/V3D.h" + +#include <boost/iterator/iterator_facade.hpp> + +namespace Mantid { +namespace Geometry { + +/** HKLGenerator + + HKLGenerator is a pseudo-container that helps in generating actual + containers with V3D-objects representing Miller indices (HKL). + + It's a common task to generate lists of Miller indices. The simplest + way of doing that is to simply create a nested loop structure that + goes through all combinations of H, K, L within some limits and then + put them into a container (vector, list, set, ...): + + for(int h = -hmin; h <= hmax; ++h) { + for(int k = -kmin; k <= kmax; ++k) { + for(int l = -lmin; l <= lmax; ++l) { + hkls.push_back(V3D(h, k, l)); + } + } + } + + In most cases the list is not used like that, instead it's filtered + using some criteria, for example a certain range of d-spacings or + others, so the reflection needs to be checked first: + + ... + hkl = V3D(h, k, l) + if(isOk(hkl)) { + hkls.push_back(hkl); + } + ... + + Instead of explicitly stating the triple-loop, HKLGenerator provides + a shorter way for this process using a const_iterator. The first code + example then becomes this: + + HKLGenerator generator(V3D(hmin, kmin, lmin), V3D(hmax, kmax, lmax)); + for(auto hkl = generator.begin(); hkl != generator.end(); ++hkl) { + hkls.push_back(*hkl); + } + + Or even shorter, using std::copy: + + HKLGenerator generator(V3D(hmin, kmin, lmin), V3D(hmax, kmax, lmax)); + std::copy(generator.begin(), generator.end(), std::back_inserter(hkls)); + + It's also possible to use filters this way, but in a much easier fashion, + using std::copy_remove_if (pre C++11) or std::copy_if (C++11): + + // pre C++11 + std::copy_remove_if(generator.begin(), generator.end(), + std::back_inserter(hkls), isNotOk) + + // C++11 + std::copy_if(generator.begin(), generator.end(), std::back_inserter(hkls), + isOk) + + See the documentation for HKLFilter for more details on how to perform + actual filtering. + + Please be aware that the iterator increments infinitely if it passes the + specified maximimum HKL. In that case K and L remain constant while H is + incremented (until it overflows). + + @author Michael Wedel, ESS + @date 23/09/2015 + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> +*/ +class MANTID_GEOMETRY_DLL HKLGenerator { +public: + /** + * @brief The const_iterator class + * + * This class uses boost::iterator_facade to provide a const_iterator + * interface to HKLGenerator. This makes it very easy to use the standard + * library algorithms. It's a forward iterator, so it's not possible to use + * operator[] for random access of specific HKLs (this would have to be + * done on the generated container if that supports it). + * + * While the iterators can be instantiated directly, the intention is to + * use HKLGenerator::begin() and HKLGenerator::end(). + */ + class const_iterator + : public boost::iterator_facade<const_iterator, const Kernel::V3D &, + boost::forward_traversal_tag> { + public: + /// Default constructor, requirement from boost::iterator_facade + const_iterator() {} + + explicit const_iterator(const Kernel::V3D ¤t); + + explicit const_iterator(const Kernel::V3D &hklMin, + const Kernel::V3D &hklMax); + + private: + // Required for boost::iterator_facade to work + friend class boost::iterator_core_access; + + void increment(); + + /// Returns true if other is at the same position + inline bool equal(const const_iterator &other) const { + return this->m_h == other.m_h && this->m_k == other.m_k && + this->m_l == other.m_l; + } + + /// Returns a const reference to the currently pointed at HKL. + inline const Kernel::V3D &dereference() const { return m_hkl; } + + int m_h, m_k, m_l; + Kernel::V3D m_hkl; + + int m_hMin, m_hMax; + int m_kMin, m_kMax; + int m_lMin, m_lMax; + }; + + HKLGenerator(const Kernel::V3D &hklMin, const Kernel::V3D &hklMax); + HKLGenerator(const Kernel::V3D &hklMinMax); + HKLGenerator(int hMinMax, int kMinMax, int lMinMax); + HKLGenerator(const UnitCell &unitCell, double dMin); + + virtual ~HKLGenerator() {} + + /// Returns the number of HKLs to be generated. + inline size_t size() const { return m_size; } + + /// Returns an iterator to the beginning of the sequence. + inline const const_iterator &begin() const { return m_begin; } + + /// Returns an iterator which "points at" one element past the end. + inline const const_iterator &end() const { return m_end; } + +private: + size_t getSize(const Kernel::V3D &min, const Kernel::V3D &max) const; + + const_iterator getBeginIterator() const; + const_iterator getEndIterator() const; + Kernel::V3D getEndHKL() const; + + Kernel::V3D m_hklMin; + Kernel::V3D m_hklMax; + + size_t m_size; + + const_iterator m_begin; + const_iterator m_end; +}; + +} // namespace Geometry +} // namespace Mantid + +#endif /* MANTID_GEOMETRY_HKLGENERATOR_H_ */ diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/ReflectionGenerator.h b/Framework/Geometry/inc/MantidGeometry/Crystal/ReflectionGenerator.h new file mode 100644 index 0000000000000000000000000000000000000000..439fcfa162daeee23ff22469fd2f258932b7ad74 --- /dev/null +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/ReflectionGenerator.h @@ -0,0 +1,144 @@ +#ifndef MANTID_GEOMETRY_REFLECTIONGENERATOR_H_ +#define MANTID_GEOMETRY_REFLECTIONGENERATOR_H_ + +#include "MantidGeometry/DllConfig.h" +#include "MantidGeometry/Crystal/CrystalStructure.h" +#include "MantidGeometry/Crystal/StructureFactorCalculator.h" +#include "MantidGeometry/Crystal/HKLFilter.h" + +namespace Mantid { +namespace Geometry { + +enum struct ReflectionConditionFilter { + None, + Centering, + SpaceGroup, + StructureFactor +}; + +/** ReflectionGenerator + + ReflectionGenerator is a class that provides the means to perform some + common tasks involving generation of reflections. While the combination + of HKLGenerator and HKLFilter is very flexible, very often a limited set + of operations has to be performed repeatedly, involving the crystal + structure. + + ReflectionGenerator is constructed from a CrystalStructure object, which + is then stored internally. Additionally, a default filter for reflection + conditions can be set, which is applied for HKL-generation in addition + to a DRangeFilter. For more flexibility, a method is provided that accepts + an HKLFilter as additional argument, this filter is then added to the + DRangeFilter. + + This way it's very simple to obtain for example a list of unique reflections + for a given crystal structure: + + CrystalStructure structure("5.43 5.43 5.43", + "F d -3 m", "Si 0 0 0 1.0 0.05"); + ReflectionGenerator generator(structure); + + // Get all unique HKLs between 0.5 and 5.0 Angstrom + std::vector<V3D> hkls = generator.getUniqueHKLs(0.5, 5.0); + + Additionally there are methods to obtain structure factors and d-values + for a given list of HKLs. + + The generated reflection lists can be filtered by several criteria to remove + reflections that are not allowed. In the default case, the reflection + conditions of the space group are used. Alternatively, the reflections can be + filtered according to the centering of the lattice or, with some more + computational effort, by their structure factors. + + The default filter method can be supplied to the constructor in the form of + the enum ReflectionConditionFilter. Furthermore, it's possible to provide + self-defined filters in overloaded versions of the HKL generation methods. + These can also be generated by using a small helper method that creates + the filters and populates them with values associated to the stored crystal + structure: + + + ReflectionGenerator generator( + CrystalStructure("5.43 5.43 5.43", + "F d -3 m", "Si 0 0 0 1.0 0.05")); + auto filter = generator.getReflectionConditionFilter( + ReflectionConditionFilter::StructureFactor); + + The filter can then be used in connection with HKLGenerator or the HKL- + generation methods of ReflectionGenerator. + + An example where structure factors are required can be found in a very + common crystal structure, the structure of Silicon. Silicon crystallizes + in the space group Fd-3m (No. 227) and the asymmetric unit consists of + one Si-atom at the position (1/8, 1/8, 1/8) (for origin choice 2, with the + inversion at the origin). Looking up this space group in the International + Tables for Crystallography A reveals that placing a scatterer at this + position introduces a new reflection condition for general reflections hkl: + + h = 2n + 1 (reflections with odd h are allowed) + or h + k + l = 4n + + This means that for example the reflection family {2 2 2} is not allowed, + even though the F-centering would allow it. Using the structure factor + calculation filter solves this problem. + + @author Michael Wedel, ESS + @date 30/09/2015 + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> +*/ +class MANTID_GEOMETRY_DLL ReflectionGenerator { +public: + ReflectionGenerator(const CrystalStructure &crystalStructure, + ReflectionConditionFilter defaultFilter = + ReflectionConditionFilter::SpaceGroup); + ~ReflectionGenerator() {} + + const CrystalStructure &getCrystalStructure() const; + + HKLFilter_const_sptr getDRangeFilter(double dMin, double dMax) const; + HKLFilter_const_sptr + getReflectionConditionFilter(ReflectionConditionFilter filter); + + std::vector<Kernel::V3D> getHKLs(double dMin, double dMax) const; + std::vector<Kernel::V3D> + getHKLs(double dMin, double dMax, + HKLFilter_const_sptr reflectionConditionFilter) const; + + std::vector<Kernel::V3D> getUniqueHKLs(double dMin, double dMax) const; + std::vector<Kernel::V3D> + getUniqueHKLs(double dMin, double dMax, + HKLFilter_const_sptr reflectionConditionFilter) const; + + std::vector<double> getDValues(const std::vector<Kernel::V3D> &hkls) const; + std::vector<double> getFsSquared(const std::vector<Kernel::V3D> &hkls) const; + +private: + CrystalStructure m_crystalStructure; + StructureFactorCalculator_sptr m_sfCalculator; + HKLFilter_const_sptr m_defaultHKLFilter; +}; + +} // namespace Geometry +} // namespace Mantid + +#endif /* MANTID_GEOMETRY_REFLECTIONGENERATOR_H_ */ diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/StructureFactorCalculator.h b/Framework/Geometry/inc/MantidGeometry/Crystal/StructureFactorCalculator.h new file mode 100644 index 0000000000000000000000000000000000000000..04b51ff0b49f016817b43884dd822c9c59e09041 --- /dev/null +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/StructureFactorCalculator.h @@ -0,0 +1,82 @@ +#ifndef MANTID_GEOMETRY_STRUCTUREFACTORCALCULATOR_H_ +#define MANTID_GEOMETRY_STRUCTUREFACTORCALCULATOR_H_ + +#include "MantidGeometry/DllConfig.h" +#include "MantidGeometry/Crystal/CrystalStructure.h" + +namespace Mantid { +namespace Geometry { + +/** StructureFactorCalculator + + This is a base class for concrete structure factor calculators. It is used + to logically separate this calculation from CrystalStructure so that different + methods of calculation can be used. For actual implementations please consult + the available sub-classes. + + @author Michael Wedel, ESS + @date 05/09/2015 + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> +*/ +class MANTID_GEOMETRY_DLL StructureFactorCalculator { +public: + StructureFactorCalculator(); + virtual ~StructureFactorCalculator() {} + + void setCrystalStructure(const CrystalStructure &crystalStructure); + + /// In implementations this method should return the structure factor for the + /// specified HKL. + virtual StructureFactor getF(const Kernel::V3D &hkl) const = 0; + virtual double getFSquared(const Kernel::V3D &hkl) const; + + virtual std::vector<StructureFactor> + getFs(const std::vector<Kernel::V3D> &hkls) const; + virtual std::vector<double> + getFsSquared(const std::vector<Kernel::V3D> &hkls) const; + +protected: + virtual void + crystalStructureSetHook(const CrystalStructure &crystalStructure); +}; + +typedef boost::shared_ptr<StructureFactorCalculator> + StructureFactorCalculator_sptr; + +namespace StructureFactorCalculatorFactory { +/// Small templated factory function that creates the desired calculator +/// and initializes it by setting the crystal structure. +template <typename T> +StructureFactorCalculator_sptr +create(const CrystalStructure &crystalStructure) { + boost::shared_ptr<T> calculator = boost::make_shared<T>(); + calculator->setCrystalStructure(crystalStructure); + + return calculator; +} +} + +} // namespace Geometry +} // namespace Mantid + +#endif /* MANTID_GEOMETRY_STRUCTUREFACTORCALCULATOR_H_ */ diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/StructureFactorCalculatorSummation.h b/Framework/Geometry/inc/MantidGeometry/Crystal/StructureFactorCalculatorSummation.h new file mode 100644 index 0000000000000000000000000000000000000000..5cfed836be0c671fee66916460800dcfc735911f --- /dev/null +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/StructureFactorCalculatorSummation.h @@ -0,0 +1,64 @@ +#ifndef MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORSUMMATION_H_ +#define MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORSUMMATION_H_ + +#include "MantidGeometry/DllConfig.h" +#include "MantidGeometry/Crystal/StructureFactorCalculator.h" + +namespace Mantid { +namespace Geometry { + +/** StructureFactorCalculatorSummation + + This implementation of StructureFactorCalculator uses the summation method + provided by BraggScatterer and its sub-classes. It obtains all scatterers in + the unit cell by combining the space group and the scatterers located in the + asymmetric unit (both taken from CrystalStructure) and stores them. + + @author Michael Wedel, ESS + @date 05/09/2015 + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> +*/ +class MANTID_GEOMETRY_DLL StructureFactorCalculatorSummation + : public StructureFactorCalculator { +public: + StructureFactorCalculatorSummation(); + virtual ~StructureFactorCalculatorSummation() {} + + StructureFactor getF(const Kernel::V3D &hkl) const; + +protected: + void crystalStructureSetHook(const CrystalStructure &crystalStructure); + + void updateUnitCellScatterers(const CrystalStructure &crystalStructure); + std::string getV3DasString(const Kernel::V3D &point) const; + + CompositeBraggScatterer_sptr m_unitCellScatterers; +}; + +typedef boost::shared_ptr<StructureFactorCalculatorSummation> + StructureFactorSummation_sptr; + +} // namespace Geometry +} // namespace Mantid + +#endif /* MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORSUMMATION_H_ */ diff --git a/Framework/Geometry/src/Crystal/BasicHKLFilters.cpp b/Framework/Geometry/src/Crystal/BasicHKLFilters.cpp new file mode 100644 index 0000000000000000000000000000000000000000..9dfec3038a8b7d458404e4ea5490392e80bc2182 --- /dev/null +++ b/Framework/Geometry/src/Crystal/BasicHKLFilters.cpp @@ -0,0 +1,120 @@ +#include "MantidGeometry/Crystal/BasicHKLFilters.h" + +namespace Mantid { +namespace Geometry { + +/// Constructor, dMax is set to the largest lattice parameter. +HKLFilterDRange::HKLFilterDRange(const UnitCell &cell, double dMin) + : m_cell(cell), m_dmin(dMin) { + m_dmax = std::max(m_cell.a(), std::max(m_cell.b(), m_cell.c())); + + checkProperDRangeValues(); +} + +/// Constructor with explicit dMax. +HKLFilterDRange::HKLFilterDRange(const UnitCell &cell, double dMin, double dMax) + : m_cell(cell), m_dmin(dMin), m_dmax(dMax) { + checkProperDRangeValues(); +} + +/// Returns a description containing the parameters of the filter. +std::string HKLFilterDRange::getDescription() const { + std::ostringstream strm; + strm << "(" << m_dmin << " <= d <= " << m_dmax << ")"; + + return strm.str(); +} + +/// Returns true if the d-value of the HKL is within the specified range. +bool HKLFilterDRange::isAllowed(const Kernel::V3D &hkl) const { + double d = m_cell.d(hkl); + + return d >= m_dmin && d <= m_dmax; +} + +/// Throws exception if m_dMin or m_dMax is <= 0 or if m_dMax < m_dMin. +void HKLFilterDRange::checkProperDRangeValues() { + if (m_dmin <= 0.0) { + throw std::range_error("dMin cannot be <= 0."); + } + + if (m_dmax <= 0.0) { + throw std::range_error("dMax cannot be <= 0."); + } + + if (m_dmax <= m_dmin) { + throw std::range_error("dMax cannot be smaller than dMin."); + } +} + +/// Constructor, throws exception if the supplied pointer is invalid. +HKLFilterSpaceGroup::HKLFilterSpaceGroup( + const SpaceGroup_const_sptr &spaceGroup) + : m_spaceGroup(spaceGroup) { + if (!m_spaceGroup) { + throw std::runtime_error( + "Cannot construct HKLFilterSpaceGroup from null space group."); + } +} + +/// Returns a description of the filter that contains the space group symbol. +std::string HKLFilterSpaceGroup::getDescription() const { + std::ostringstream strm; + strm << "(Space group: " << m_spaceGroup->hmSymbol() << ")"; + + return strm.str(); +} + +/// Returns true if the reflection is allowed by the space group reflection +/// conditions. +bool HKLFilterSpaceGroup::isAllowed(const Kernel::V3D &hkl) const { + return m_spaceGroup->isAllowedReflection(hkl); +} + +/// Constructor, throws exception if the calculator pointer is invalid. +HKLFilterStructureFactor::HKLFilterStructureFactor( + const StructureFactorCalculator_sptr &calculator, double fSquaredMin) + : m_calculator(calculator), m_fSquaredMin(fSquaredMin) { + if (!m_calculator) { + throw std::runtime_error( + "Cannot construct HKLFilterStructureFactor from null calculator."); + } +} + +/// Returns a description for the filter that contains the minimum F^2. +std::string HKLFilterStructureFactor::getDescription() const { + std::ostringstream strm; + strm << "(F^2 > " << m_fSquaredMin << ")"; + + return strm.str(); +} + +/// Returns true if F^2(hkl) is larger than the stored minimum. +bool HKLFilterStructureFactor::isAllowed(const Kernel::V3D &hkl) const { + return m_calculator->getFSquared(hkl) > m_fSquaredMin; +} + +/// Constructor, throws exception if pointer is null. +HKLFilterCentering::HKLFilterCentering( + const ReflectionCondition_sptr ¢ering) + : m_centering(centering) { + if (!m_centering) { + throw std::runtime_error( + "Cannot construct HKLFilterCentering from null centering."); + } +} + +/// Returns a description with the centering symbol. +std::string HKLFilterCentering::getDescription() const { + return "(Centering: " + m_centering->getSymbol() + ")"; +} + +/// Returns true if the HKL is allowed according to the lattice centering. +bool HKLFilterCentering::isAllowed(const Kernel::V3D &hkl) const { + return m_centering->isAllowed(static_cast<int>(hkl.X()), + static_cast<int>(hkl.Y()), + static_cast<int>(hkl.Z())); +} + +} // namespace Geometry +} // namespace Mantid diff --git a/Framework/Geometry/src/Crystal/HKLFilter.cpp b/Framework/Geometry/src/Crystal/HKLFilter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ed218cb7fdb2f3fce2bdeaa204efc113e8479009 --- /dev/null +++ b/Framework/Geometry/src/Crystal/HKLFilter.cpp @@ -0,0 +1,117 @@ +#include "MantidGeometry/Crystal/HKLFilter.h" +#include <boost/make_shared.hpp> +#include <stdexcept> + +namespace Mantid { +namespace Geometry { + +/** + * Returns a function object that wraps HKLFilter::isAllowed + * + * This method uses std::bind to create a function object that + * represents the HKLFilter::isAllowed() method. This way it's + * possible to pass the function to STL-algorithms easily (see + * class documentation). + * + * @return Function object with filter function for V3D.s + */ +std::function<bool(const Kernel::V3D &)> HKLFilter::fn() const { + return std::bind(&HKLFilter::isAllowed, this, std::placeholders::_1); +} + +/// Stores the supplied filter, throws exception if filter is null. +HKLFilterUnaryLogicOperation::HKLFilterUnaryLogicOperation( + const HKLFilter_const_sptr &filter) + : m_operand(filter) { + if (!m_operand) { + throw std::runtime_error( + "Cannot create HKLFilterUnaryLogicOperation from null operand."); + } +} + +/// Returns a description of the HKLFilterNot. +std::string HKLFilterNot::getDescription() const { + return "!" + m_operand->getDescription(); +} + +/// Returns true if the wrapped filter returns false and false otherwise. +bool HKLFilterNot::isAllowed(const Kernel::V3D &hkl) const { + return !(m_operand->isAllowed(hkl)); +} + +/// Stores the left-hand and right-hand side operators, throws exception if +/// either is null. +HKLFilterBinaryLogicOperation::HKLFilterBinaryLogicOperation( + const HKLFilter_const_sptr &lhs, const HKLFilter_const_sptr &rhs) + : m_lhs(lhs), m_rhs(rhs) { + if (!m_lhs || !m_rhs) { + throw std::runtime_error("Cannot construct HKLFilterBinaryLogicOperation " + "with one or more null-operands."); + } +} + +/// Returns a description of the HKLFilterAnd. +std::string HKLFilterAnd::getDescription() const { + return "(" + m_lhs->getDescription() + " & " + m_rhs->getDescription() + ")"; +} + +/// Returns true if both wrapped filters return true. +bool HKLFilterAnd::isAllowed(const Kernel::V3D &hkl) const { + return m_lhs->isAllowed(hkl) && m_rhs->isAllowed(hkl); +} + +/// Returns a description of the HKLFilterOr. +std::string HKLFilterOr::getDescription() const { + return "(" + m_lhs->getDescription() + " | " + m_rhs->getDescription() + ")"; +} + +/// Returns true if either of the wrapped filters returns true. +bool HKLFilterOr::isAllowed(const Kernel::V3D &hkl) const { + return m_lhs->isAllowed(hkl) || m_rhs->isAllowed(hkl); +} + +/** + * Constructs an HKLFilterNot from the operand + * + * This function makes it easy to construct HKLFilterNot by using the + * operator ~, which inverts the wrapped filter. + * + * @param filter :: HKLFilter to invert. + * @return HKLFilterNot with the wrapped filter. + */ +const HKLFilter_const_sptr operator~(const HKLFilter_const_sptr &filter) { + return boost::make_shared<const HKLFilterNot>(filter); +} + +/** + * Constructs an HKLFilterAnd from the operands + * + * This function makes it easy to construct HKLFilterAnd by using the + * operator & directly on HKLFilter_const_sptr. + * + * @param lhs :: Left-hand side HKLFilter operand. + * @param rhs :: Right-hand side HKLFilter operand. + * @return HKLFilterAnd with two wrapped filters. + */ +const HKLFilter_const_sptr operator&(const HKLFilter_const_sptr &lhs, + const HKLFilter_const_sptr &rhs) { + return boost::make_shared<const HKLFilterAnd>(lhs, rhs); +} + +/** + * Constructs an HKLFilterOr from the operands + * + * This function makes it easy to construct HKLFilterOr by using the + * operator | directly on HKLFilter_const_sptr. + * + * @param lhs :: Left-hand side HKLFilter operand. + * @param rhs :: Right-hand side HKLFilter operand. + * @return HKLFilterOr with two wrapped filters. + */ +const HKLFilter_const_sptr operator|(const HKLFilter_const_sptr &lhs, + const HKLFilter_const_sptr &rhs) { + return boost::make_shared<HKLFilterOr>(lhs, rhs); +} + +} // namespace Geometry +} // namespace Mantid diff --git a/Framework/Geometry/src/Crystal/HKLGenerator.cpp b/Framework/Geometry/src/Crystal/HKLGenerator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0156cd97045b0a347af482893774907c011935c0 --- /dev/null +++ b/Framework/Geometry/src/Crystal/HKLGenerator.cpp @@ -0,0 +1,96 @@ +#include "MantidGeometry/Crystal/HKLGenerator.h" + +namespace Mantid { +namespace Geometry { + +using namespace Kernel; + +/// Constructs a generator that creates all indices from hklMin to hklMax. +HKLGenerator::HKLGenerator(const Kernel::V3D &hklMin, const Kernel::V3D &hklMax) + : m_hklMin(hklMin), m_hklMax(hklMax), m_size(getSize(m_hklMin, m_hklMax)), + m_begin(getBeginIterator()), m_end(getEndIterator()) {} + +/// Constructs a generator that creates all indices from -hklMinMax to +/// hklMinMax. +HKLGenerator::HKLGenerator(const Kernel::V3D &hklMinMax) + : m_hklMin(hklMinMax * -1), m_hklMax(hklMinMax), + m_size(getSize(m_hklMin, m_hklMax)), m_begin(getBeginIterator()), + m_end(getEndIterator()) {} + +/// Constructs a generator that creates all indices from -h,-k,-l to h,k,l. +HKLGenerator::HKLGenerator(int hMinMax, int kMinMax, int lMinMax) { + m_hklMax = V3D(hMinMax, kMinMax, lMinMax); + m_hklMin = m_hklMax * -1; + m_size = getSize(m_hklMin, m_hklMax); + + m_begin = getBeginIterator(); + m_end = getEndIterator(); +} + +/// Constructs a generator that creates all indices for the given cell up to +/// dMin. +HKLGenerator::HKLGenerator(const UnitCell &unitCell, double dMin) { + m_hklMax = V3D(unitCell.a() / dMin, unitCell.b() / dMin, unitCell.c() / dMin); + m_hklMin = m_hklMax * -1; + m_size = getSize(m_hklMin, m_hklMax); + + m_begin = getBeginIterator(); + m_end = getEndIterator(); +} + +/// Returns the number of indices between min and max. +size_t HKLGenerator::getSize(const V3D &min, const V3D &max) const { + V3D diff = (max - min) + V3D(1, 1, 1); + return static_cast<size_t>(diff.X() * diff.Y() * diff.Z()); +} + +/// Constructs an iterator that points to the beginning of the sequence. +HKLGenerator::const_iterator HKLGenerator::getBeginIterator() const { + return HKLGenerator::const_iterator(m_hklMin, m_hklMax); +} + +/// Constructs an iterator that points to an HKL one past the maximum. +HKLGenerator::const_iterator HKLGenerator::getEndIterator() const { + return HKLGenerator::const_iterator(getEndHKL()); +} + +/// Returns the HKL "one past the maximum". +V3D HKLGenerator::getEndHKL() const { + return V3D(m_hklMax.X() + 1, m_hklMin.Y(), m_hklMin.Z()); +} + +/// Return an iterator with min = max = current. +HKLGenerator::const_iterator::const_iterator(const V3D ¤t) + : m_h(static_cast<int>(current.X())), m_k(static_cast<int>(current.Y())), + m_l(static_cast<int>(current.Z())), m_hkl(current), m_hMin(m_h), + m_hMax(m_h), m_kMin(m_k), m_kMax(m_k), m_lMin(m_l), m_lMax(m_l) {} + +/// Return an iterator that can move from min to max, with current = min +HKLGenerator::const_iterator::const_iterator(const V3D &hklMin, + const V3D &hklMax) + : m_h(static_cast<int>(hklMin.X())), m_k(static_cast<int>(hklMin.Y())), + m_l(static_cast<int>(hklMin.Z())), m_hkl(hklMin), m_hMin(m_h), + m_hMax(static_cast<int>(hklMax.X())), m_kMin(m_k), + m_kMax(static_cast<int>(hklMax.Y())), m_lMin(m_l), + m_lMax(static_cast<int>(hklMax.Z())) {} + +/// Increments HKL, l moves fastest, h moves slowest, wrapping around at the max +/// for each index. +void HKLGenerator::const_iterator::increment() { + ++m_l; + + if (m_l > m_lMax) { + m_l = m_lMin; + + ++m_k; + if (m_k > m_kMax) { + m_k = m_kMin; + ++m_h; + } + } + + m_hkl = V3D(m_h, m_k, m_l); +} + +} // namespace Geometry +} // namespace Mantid diff --git a/Framework/Geometry/src/Crystal/ReflectionGenerator.cpp b/Framework/Geometry/src/Crystal/ReflectionGenerator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..59fd597ef216ad55f5d601c3a8f688c5e443738d --- /dev/null +++ b/Framework/Geometry/src/Crystal/ReflectionGenerator.cpp @@ -0,0 +1,146 @@ +#include "MantidGeometry/Crystal/ReflectionGenerator.h" +#include "MantidGeometry/Crystal/BasicHKLFilters.h" +#include "MantidGeometry/Crystal/StructureFactorCalculatorSummation.h" +#include "MantidGeometry/Crystal/HKLGenerator.h" + +namespace Mantid { +namespace Geometry { + +using namespace Kernel; + +/// Small helper functor to calculate d-Values from a unit cell. +class LatticeSpacingCalculator { +public: + LatticeSpacingCalculator(const UnitCell &cell) : m_cell(cell) {} + + double operator()(const V3D &hkl) { return m_cell.d(hkl); } + +private: + UnitCell m_cell; +}; + +/// Constructor +ReflectionGenerator::ReflectionGenerator( + const CrystalStructure &crystalStructure, + ReflectionConditionFilter defaultFilter) + : m_crystalStructure(crystalStructure), + m_sfCalculator(StructureFactorCalculatorFactory::create< + StructureFactorCalculatorSummation>(m_crystalStructure)), + m_defaultHKLFilter(getReflectionConditionFilter(defaultFilter)) {} + +/// Returns the internally stored crystal structure +const CrystalStructure &ReflectionGenerator::getCrystalStructure() const { + return m_crystalStructure; +} + +/// Returns a DRangeFilter from the supplied d-limits and the internally stored +/// cell. +HKLFilter_const_sptr ReflectionGenerator::getDRangeFilter(double dMin, + double dMax) const { + return boost::make_shared<const HKLFilterDRange>(m_crystalStructure.cell(), + dMin, dMax); +} + +/// Returns a reflection condition HKLFilter based on the supplied enum. +HKLFilter_const_sptr ReflectionGenerator::getReflectionConditionFilter( + ReflectionConditionFilter filter) { + switch (filter) { + case ReflectionConditionFilter::Centering: + return boost::make_shared<const HKLFilterCentering>( + m_crystalStructure.centering()); + break; + case ReflectionConditionFilter::SpaceGroup: + return boost::make_shared<const HKLFilterSpaceGroup>( + m_crystalStructure.spaceGroup()); + break; + case ReflectionConditionFilter::StructureFactor: + return boost::make_shared<const HKLFilterStructureFactor>(m_sfCalculator); + default: + return HKLFilter_const_sptr(); + } +} + +/// Returns a list of HKLs within the specified d-limits using the default +/// reflection condition filter. +std::vector<V3D> ReflectionGenerator::getHKLs(double dMin, double dMax) const { + return getHKLs(dMin, dMax, m_defaultHKLFilter); +} + +/// Returns a list of HKLs within the specified d-limits using the specified +/// filter. If the pointer is null, it's ignored. +std::vector<Kernel::V3D> ReflectionGenerator::getHKLs( + double dMin, double dMax, + HKLFilter_const_sptr reflectionConditionFilter) const { + HKLGenerator generator(m_crystalStructure.cell(), dMin); + + HKLFilter_const_sptr filter = getDRangeFilter(dMin, dMax); + if (reflectionConditionFilter) { + filter = filter & reflectionConditionFilter; + } + + std::vector<V3D> hkls; + hkls.reserve(generator.size()); + + std::remove_copy_if(generator.begin(), generator.end(), + std::back_inserter(hkls), (~filter)->fn()); + return hkls; +} + +/// Returns a list of symetrically independent HKLs within the specified +/// d-limits using the default reflection condition filter. +std::vector<V3D> ReflectionGenerator::getUniqueHKLs(double dMin, + double dMax) const { + return getUniqueHKLs(dMin, dMax, m_defaultHKLFilter); +} + +/// Returns a list of symetrically independent HKLs within the specified +/// d-limits using the specified reflection condition filter. +std::vector<V3D> ReflectionGenerator::getUniqueHKLs( + double dMin, double dMax, + HKLFilter_const_sptr reflectionConditionFilter) const { + HKLGenerator generator(m_crystalStructure.cell(), dMin); + + HKLFilter_const_sptr filter = getDRangeFilter(dMin, dMax); + if (reflectionConditionFilter) { + filter = filter & reflectionConditionFilter; + } + + std::vector<V3D> hkls; + hkls.reserve(generator.size()); + + PointGroup_sptr pg = m_crystalStructure.spaceGroup()->getPointGroup(); + + for (auto hkl = generator.begin(); hkl != generator.end(); ++hkl) { + if (filter->isAllowed(*hkl)) { + hkls.push_back(pg->getReflectionFamily(*hkl)); + } + } + + std::sort(hkls.begin(), hkls.end()); + hkls.erase(std::unique(hkls.begin(), hkls.end()), hkls.end()); + + return hkls; +} + +/// Returns a list of d-values that correspond to the supplied hkl list, using +/// the unit cell of the stored crystal structure. +std::vector<double> +ReflectionGenerator::getDValues(const std::vector<V3D> &hkls) const { + std::vector<double> dValues; + dValues.reserve(hkls.size()); + + std::transform(hkls.begin(), hkls.end(), std::back_inserter(dValues), + LatticeSpacingCalculator(m_crystalStructure.cell())); + + return dValues; +} + +/// Returns a list of squared structure factor amplitudes corresponding to the +/// supplied list of HKLs. +std::vector<double> +ReflectionGenerator::getFsSquared(const std::vector<V3D> &hkls) const { + return m_sfCalculator->getFsSquared(hkls); +} + +} // namespace Geometry +} // namespace Mantid diff --git a/Framework/Geometry/src/Crystal/StructureFactorCalculator.cpp b/Framework/Geometry/src/Crystal/StructureFactorCalculator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..df4cb44ed97825cee453bb450b2be33a9b5be755 --- /dev/null +++ b/Framework/Geometry/src/Crystal/StructureFactorCalculator.cpp @@ -0,0 +1,83 @@ +#include "MantidGeometry/Crystal/StructureFactorCalculator.h" + +#include <boost/bind.hpp> + +namespace Mantid { +namespace Geometry { + +StructureFactorCalculator::StructureFactorCalculator() {} + +/** + * Sets the crystal structure for which to calculate structure factors + * + * Additionally StructureFactorCalculator::crystalStructureSetHook() is called + * with the CrystalStructure. This function may be re-implemented + * by concrete structure factor calculators to perform other necessary actions + * that involve the crystal structure (for example extracting and storing + * members of it). + * + * @param crystalStructure :: Crystal structure for calculations. + */ +void StructureFactorCalculator::setCrystalStructure( + const CrystalStructure &crystalStructure) { + crystalStructureSetHook(crystalStructure); +} + +/// Returns F^2 for the given HKL, calling StructureFactorCalculator::getF(). +double StructureFactorCalculator::getFSquared(const Kernel::V3D &hkl) const { + StructureFactor sf = getF(hkl); + + return sf.real() * sf.real() + sf.imag() * sf.imag(); +} + +/** + * Returns structure factors for each HKL in the container + * + * The default implementation uses StructureFactorCalculator::getF() to get + * the structure factor for each HKL. This behavior can be overriden in + * sub-classes for example to implement this in a more efficient way. + * + * @param hkls :: Vector of HKLs. + * @return :: Vector of structure factors for the given HKLs. + */ +std::vector<StructureFactor> +StructureFactorCalculator::getFs(const std::vector<Kernel::V3D> &hkls) const { + std::vector<StructureFactor> structureFactors(hkls.size()); + + std::transform(hkls.begin(), hkls.end(), structureFactors.begin(), + boost::bind(&StructureFactorCalculator::getF, this, _1)); + + return structureFactors; +} + +/** + * Returns structure factors for each HKL in the container + * + * The default implementation uses StructureFactorCalculator::getFSquared() with + * each HKL. This behavior can be overriden in + * sub-classes for example to implement this in a more efficient way. + * + * @param hkls :: Vector of HKLs. + * @return :: Vector of squared structure factors for the given HKLs. + */ +std::vector<double> StructureFactorCalculator::getFsSquared( + const std::vector<Kernel::V3D> &hkls) const { + std::vector<double> fSquareds(hkls.size()); + + std::transform( + hkls.begin(), hkls.end(), fSquareds.begin(), + boost::bind(&StructureFactorCalculator::getFSquared, this, _1)); + + return fSquareds; +} + +/// This function is called from +/// StructureFactorCalculator::setCrystalStructure() and can be overriden to +/// perform additional actions. +void StructureFactorCalculator::crystalStructureSetHook( + const CrystalStructure &crystalStructure) { + UNUSED_ARG(crystalStructure) +} + +} // namespace Geometry +} // namespace Mantid diff --git a/Framework/Geometry/src/Crystal/StructureFactorCalculatorSummation.cpp b/Framework/Geometry/src/Crystal/StructureFactorCalculatorSummation.cpp new file mode 100644 index 0000000000000000000000000000000000000000..eabc22952cbd37217570255a766d4f023f170d83 --- /dev/null +++ b/Framework/Geometry/src/Crystal/StructureFactorCalculatorSummation.cpp @@ -0,0 +1,82 @@ +#include "MantidGeometry/Crystal/StructureFactorCalculatorSummation.h" +#include "MantidGeometry/Crystal/BraggScattererInCrystalStructure.h" + +#include <iomanip> + +namespace Mantid { +namespace Geometry { + +using namespace Kernel; + +StructureFactorCalculatorSummation::StructureFactorCalculatorSummation() + : StructureFactorCalculator(), + m_unitCellScatterers(CompositeBraggScatterer::create()) {} + +/// Returns the structure factor obtained from the stored scatterers. +StructureFactor +StructureFactorCalculatorSummation::getF(const Kernel::V3D &hkl) const { + return m_unitCellScatterers->calculateStructureFactor(hkl); +} + +/// Calls updateUnitCellScatterers() to rebuild the complete list of scatterers. +void StructureFactorCalculatorSummation::crystalStructureSetHook( + const CrystalStructure &crystalStructure) { + updateUnitCellScatterers(crystalStructure); +} + +/** + * Rebuilds the internal list of scatterers + * + * This method extracts two items from the supplied CrystalStructure object, the + * list of scatterers in the asymmetric unit and the space group. Then it + * generates all equivalent positions and creates a complete list of scatterers. + * + * @param crystalStructure :: CrystalStructure for structure factor calculation. + */ +void StructureFactorCalculatorSummation::updateUnitCellScatterers( + const CrystalStructure &crystalStructure) { + m_unitCellScatterers->removeAllScatterers(); + + CompositeBraggScatterer_sptr scatterersInAsymmetricUnit = + crystalStructure.getScatterers(); + SpaceGroup_const_sptr spaceGroup = crystalStructure.spaceGroup(); + + if (spaceGroup) { + std::vector<BraggScatterer_sptr> braggScatterers; + braggScatterers.reserve(scatterersInAsymmetricUnit->nScatterers() * + spaceGroup->order()); + + for (size_t i = 0; i < scatterersInAsymmetricUnit->nScatterers(); ++i) { + BraggScattererInCrystalStructure_sptr current = + boost::dynamic_pointer_cast<BraggScattererInCrystalStructure>( + scatterersInAsymmetricUnit->getScatterer(i)); + + if (current) { + std::vector<V3D> positions = + spaceGroup->getEquivalentPositions(current->getPosition()); + + for (auto pos = positions.begin(); pos != positions.end(); ++pos) { + BraggScatterer_sptr clone = current->clone(); + clone->setProperty("Position", getV3DasString(*pos)); + + braggScatterers.push_back(clone); + } + } + } + + m_unitCellScatterers->setScatterers(braggScatterers); + } +} + +/// Return V3D as string without losing precision. +std::string +StructureFactorCalculatorSummation::getV3DasString(const V3D &point) const { + std::ostringstream posStream; + posStream << std::setprecision(17); + posStream << point; + + return posStream.str(); +} + +} // namespace Geometry +} // namespace Mantid diff --git a/Framework/Geometry/test/BasicHKLFiltersTest.h b/Framework/Geometry/test/BasicHKLFiltersTest.h new file mode 100644 index 0000000000000000000000000000000000000000..7295290c09a156062be0e9cd70d3c2308ef2bead --- /dev/null +++ b/Framework/Geometry/test/BasicHKLFiltersTest.h @@ -0,0 +1,146 @@ +#ifndef MANTID_GEOMETRY_BASICHKLFILTERSTEST_H_ +#define MANTID_GEOMETRY_BASICHKLFILTERSTEST_H_ + +#include <cxxtest/TestSuite.h> +#include <gtest/gtest.h> +#include <gmock/gmock.h> + +#include "MantidGeometry/Crystal/BasicHKLFilters.h" +#include "MantidGeometry/Crystal/SpaceGroupFactory.h" +#include "MantidGeometry/Crystal/StructureFactorCalculator.h" + +#include <strstream> + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; + +using ::testing::_; +using ::testing::Return; +using ::testing::Mock; + +class BasicHKLFiltersTest : public CxxTest::TestSuite { +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static BasicHKLFiltersTest *createSuite() { + return new BasicHKLFiltersTest(); + } + static void destroySuite(BasicHKLFiltersTest *suite) { delete suite; } + + void testHKLFilterDRangeConstructors() { + UnitCell cell(10., 10., 10.); + + TS_ASSERT_THROWS_NOTHING(HKLFilterDRange dFilter(cell, 1.0)); + TS_ASSERT_THROWS(HKLFilterDRange dFilter(cell, -1.0), std::range_error); + TS_ASSERT_THROWS(HKLFilterDRange dFilter(cell, 0.0), std::range_error); + + TS_ASSERT_THROWS_NOTHING(HKLFilterDRange dFilter(cell, 1.0, 2.0)); + TS_ASSERT_THROWS(HKLFilterDRange dFilter(cell, 1.0, 0.5), std::range_error); + TS_ASSERT_THROWS(HKLFilterDRange dFilter(cell, 1.0, -0.5), + std::range_error); + } + + void testHKLFilterDRangeDescription() { + UnitCell cell(10., 10., 10.); + + std::ostringstream reference; + reference << "(" << 1.0 << " <= d <= " << 10.0 << ")"; + + HKLFilterDRange dFilter(cell, 1.0); + TS_ASSERT_EQUALS(dFilter.getDescription(), reference.str()); + } + + void testHKLFilterDRangeIsAllowed() { + UnitCell cell(10., 10., 10.); + HKLFilterDRange dFilter(cell, 1.0, 9.0); + + TS_ASSERT(dFilter.isAllowed(V3D(1, 2, 3))); + + TS_ASSERT(dFilter.isAllowed(V3D(2, 0, 0))); + TS_ASSERT(!dFilter.isAllowed(V3D(1, 0, 0))); + + TS_ASSERT(dFilter.isAllowed(V3D(10, 0, 0))); + TS_ASSERT(!dFilter.isAllowed(V3D(11, 0, 0))); + } + + void testHKLFilterSpaceGroupConstructor() { + SpaceGroup_const_sptr invalid; + TS_ASSERT_THROWS(HKLFilterSpaceGroup sgFilter(invalid), std::runtime_error); + + SpaceGroup_const_sptr sg = + SpaceGroupFactory::Instance().createSpaceGroup("F d -3 m"); + TS_ASSERT_THROWS_NOTHING(HKLFilterSpaceGroup sgFilter(sg)); + } + + void testHKLFilterSpaceGroupDescription() { + SpaceGroup_const_sptr sg = + SpaceGroupFactory::Instance().createSpaceGroup("F d -3 m"); + + HKLFilterSpaceGroup sgFilter(sg); + + TS_ASSERT_EQUALS(sgFilter.getDescription(), + "(Space group: " + sg->hmSymbol() + ")"); + } + + void testHKLFilterSpaceGroupIsAllowed() { + SpaceGroup_const_sptr sg = + SpaceGroupFactory::Instance().createSpaceGroup("F d -3 m"); + + HKLFilterSpaceGroup sgFilter(sg); + + TS_ASSERT(!sgFilter.isAllowed(V3D(1, 0, 0))); + TS_ASSERT(!sgFilter.isAllowed(V3D(1, 1, 0))); + TS_ASSERT(sgFilter.isAllowed(V3D(1, 1, 1))); + + TS_ASSERT(!sgFilter.isAllowed(V3D(2, 0, 0))); + TS_ASSERT(!sgFilter.isAllowed(V3D(3, 0, 0))); + TS_ASSERT(sgFilter.isAllowed(V3D(4, 0, 0))); + } + + void testHKLFilterStructureFactorConstructor() { + StructureFactorCalculator_sptr invalid; + TS_ASSERT_THROWS(HKLFilterStructureFactor sfFilter(invalid), + std::runtime_error); + + StructureFactorCalculator_sptr mock = + boost::make_shared<MockStructureFactorCalculator>(); + TS_ASSERT_THROWS_NOTHING(HKLFilterStructureFactor sfFilter(mock)); + TS_ASSERT_THROWS_NOTHING(HKLFilterStructureFactor sfFilter(mock, 12.0)); + } + + void testHKLFilterStructureFactorDescription() { + std::ostringstream reference; + reference << "(F^2 > " << 1.0 << ")"; + + StructureFactorCalculator_sptr mock = + boost::make_shared<MockStructureFactorCalculator>(); + HKLFilterStructureFactor sfFilter(mock, 1.0); + TS_ASSERT_EQUALS(sfFilter.getDescription(), reference.str()); + } + + void testHKLFilterStructureFactorIsAllowed() { + boost::shared_ptr<MockStructureFactorCalculator> mock = + boost::make_shared<MockStructureFactorCalculator>(); + + EXPECT_CALL(*mock, getFSquared(_)) + .WillOnce(Return(2.0)) + .WillOnce(Return(0.5)) + .WillOnce(Return(1.0)); + + HKLFilterStructureFactor sfFilter(mock, 1.0); + TS_ASSERT(sfFilter.isAllowed(V3D(1, 1, 1))); + TS_ASSERT(!sfFilter.isAllowed(V3D(1, 1, 1))); + TS_ASSERT(!sfFilter.isAllowed(V3D(1, 1, 1))); + + TS_ASSERT(Mock::VerifyAndClearExpectations(mock.get())) + } + +private: + class MockStructureFactorCalculator : public StructureFactorCalculator { + public: + MOCK_CONST_METHOD1(getF, StructureFactor(const V3D &)); + MOCK_CONST_METHOD1(getFSquared, double(const V3D &)); + }; +}; + +#endif /* MANTID_GEOMETRY_BASICHKLFILTERSTEST_H_ */ diff --git a/Framework/Geometry/test/HKLFilterTest.h b/Framework/Geometry/test/HKLFilterTest.h new file mode 100644 index 0000000000000000000000000000000000000000..c0ef32dde759a8601d65ffdfd9fae6adc7e001e8 --- /dev/null +++ b/Framework/Geometry/test/HKLFilterTest.h @@ -0,0 +1,209 @@ +#ifndef MANTID_GEOMETRY_HKLFILTERTEST_H_ +#define MANTID_GEOMETRY_HKLFILTERTEST_H_ + +#include <cxxtest/TestSuite.h> +#include <gtest/gtest.h> +#include <gmock/gmock.h> + +#include "MantidGeometry/Crystal/HKLFilter.h" +#include "MantidKernel/V3D.h" + +#include <boost/make_shared.hpp> + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; + +using ::testing::_; +using ::testing::Return; +using ::testing::Mock; + +class HKLFilterTest : public CxxTest::TestSuite { +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static HKLFilterTest *createSuite() { return new HKLFilterTest(); } + static void destroySuite(HKLFilterTest *suite) { delete suite; } + + void testFn() { + MockHKLFilter filter; + EXPECT_CALL(filter, isAllowed(_)) + .Times(2) + .WillOnce(Return(true)) + .WillRepeatedly(Return(false)); + + std::function<bool(const V3D &)> f = filter.fn(); + + TS_ASSERT_EQUALS(f(V3D(1, 1, 1)), true); + TS_ASSERT_EQUALS(f(V3D(1, 1, 1)), false); + + TS_ASSERT(Mock::VerifyAndClearExpectations(&filter)); + } + + void testUnaryLogicOperation() { + HKLFilter_const_sptr filter = boost::make_shared<MockHKLFilter>(); + + TS_ASSERT_THROWS_NOTHING(MockHKLFilterUnaryLogicOperation op(filter)); + + MockHKLFilterUnaryLogicOperation op(filter); + TS_ASSERT_EQUALS(op.getOperand(), filter); + + HKLFilter_const_sptr invalid; + TS_ASSERT_THROWS(MockHKLFilterUnaryLogicOperation op(invalid), + std::runtime_error); + } + + void testBinaryLogicOperation() { + HKLFilter_const_sptr lhs = boost::make_shared<MockHKLFilter>(); + HKLFilter_const_sptr rhs = boost::make_shared<MockHKLFilter>(); + + TS_ASSERT_THROWS_NOTHING(MockHKLFilterBinaryLogicOperation op(lhs, rhs)); + + MockHKLFilterBinaryLogicOperation op(lhs, rhs); + TS_ASSERT_EQUALS(op.getLHS(), lhs); + TS_ASSERT_EQUALS(op.getRHS(), rhs); + + HKLFilter_const_sptr invalid; + TS_ASSERT_THROWS(MockHKLFilterBinaryLogicOperation op(invalid, rhs), + std::runtime_error); + TS_ASSERT_THROWS(MockHKLFilterBinaryLogicOperation op(lhs, invalid), + std::runtime_error); + TS_ASSERT_THROWS(MockHKLFilterBinaryLogicOperation op(invalid, invalid), + std::runtime_error); + } + + void testHKLFilterNot() { + boost::shared_ptr<const MockHKLFilter> filter(new MockHKLFilter); + + EXPECT_CALL(*filter, isAllowed(_)) + .WillOnce(Return(true)) + .WillOnce(Return(false)); + + HKLFilterNot notFilter(filter); + + TS_ASSERT_EQUALS(notFilter.isAllowed(V3D(1, 1, 1)), false); + TS_ASSERT_EQUALS(notFilter.isAllowed(V3D(1, 1, 1)), true); + + TS_ASSERT(Mock::VerifyAndClearExpectations( + boost::const_pointer_cast<MockHKLFilter>(filter).get())); + } + + void testHKLFilterNotOperator() { + HKLFilter_const_sptr filter = boost::make_shared<MockHKLFilter>(); + + HKLFilter_const_sptr notFilter = ~filter; + + boost::shared_ptr<const HKLFilterNot> notFilterCasted = + boost::dynamic_pointer_cast<const HKLFilterNot>(notFilter); + TS_ASSERT(notFilterCasted); + + TS_ASSERT_EQUALS(notFilterCasted->getOperand(), filter); + } + + void testHKLFilterAnd() { + boost::shared_ptr<const MockHKLFilter> lhs(new MockHKLFilter); + EXPECT_CALL(*lhs, isAllowed(_)) + .WillOnce(Return(true)) + .WillOnce(Return(false)) + .WillOnce(Return(true)); + + boost::shared_ptr<const MockHKLFilter> rhs(new MockHKLFilter); + EXPECT_CALL(*rhs, isAllowed(_)) + .WillOnce(Return(true)) + .WillOnce(Return(false)); + + HKLFilterAnd andFilter(lhs, rhs); + + TS_ASSERT_EQUALS(andFilter.isAllowed(V3D(1, 1, 1)), true); + TS_ASSERT_EQUALS(andFilter.isAllowed(V3D(1, 1, 1)), false); + TS_ASSERT_EQUALS(andFilter.isAllowed(V3D(1, 1, 1)), false); + + TS_ASSERT(Mock::VerifyAndClearExpectations( + boost::const_pointer_cast<MockHKLFilter>(lhs).get())); + TS_ASSERT(Mock::VerifyAndClearExpectations( + boost::const_pointer_cast<MockHKLFilter>(rhs).get())); + } + + void testHKLFilterAndOperator() { + HKLFilter_const_sptr lhs = boost::make_shared<MockHKLFilter>(); + HKLFilter_const_sptr rhs = boost::make_shared<MockHKLFilter>(); + + HKLFilter_const_sptr andFilter = lhs & rhs; + + boost::shared_ptr<const HKLFilterAnd> andFilterCasted = + boost::dynamic_pointer_cast<const HKLFilterAnd>(andFilter); + + TS_ASSERT(andFilterCasted); + + TS_ASSERT_EQUALS(andFilterCasted->getLHS(), lhs); + TS_ASSERT_EQUALS(andFilterCasted->getRHS(), rhs); + } + + void testHKLFilterOr() { + boost::shared_ptr<const MockHKLFilter> lhs(new MockHKLFilter); + EXPECT_CALL(*lhs, isAllowed(_)) + .WillOnce(Return(true)) + .WillOnce(Return(false)) + .WillOnce(Return(true)) + .WillOnce(Return(false)); + + boost::shared_ptr<const MockHKLFilter> rhs(new MockHKLFilter); + EXPECT_CALL(*rhs, isAllowed(_)) + .WillOnce(Return(false)) + .WillOnce(Return(true)); + + HKLFilterOr orFilter(lhs, rhs); + + TS_ASSERT_EQUALS(orFilter.isAllowed(V3D(1, 1, 1)), true); + TS_ASSERT_EQUALS(orFilter.isAllowed(V3D(1, 1, 1)), false); + TS_ASSERT_EQUALS(orFilter.isAllowed(V3D(1, 1, 1)), true); + TS_ASSERT_EQUALS(orFilter.isAllowed(V3D(1, 1, 1)), true); + + TS_ASSERT(Mock::VerifyAndClearExpectations( + boost::const_pointer_cast<MockHKLFilter>(lhs).get())); + TS_ASSERT(Mock::VerifyAndClearExpectations( + boost::const_pointer_cast<MockHKLFilter>(rhs).get())); + } + + void testHKLFilterOrOperator() { + HKLFilter_const_sptr lhs = boost::make_shared<MockHKLFilter>(); + HKLFilter_const_sptr rhs = boost::make_shared<MockHKLFilter>(); + + HKLFilter_const_sptr orFilter = lhs | rhs; + + boost::shared_ptr<const HKLFilterOr> orFilterCasted = + boost::dynamic_pointer_cast<const HKLFilterOr>(orFilter); + + TS_ASSERT(orFilterCasted); + + TS_ASSERT_EQUALS(orFilterCasted->getLHS(), lhs); + TS_ASSERT_EQUALS(orFilterCasted->getRHS(), rhs); + } + +private: + class MockHKLFilter : public HKLFilter { + public: + MOCK_CONST_METHOD0(getDescription, std::string()); + MOCK_CONST_METHOD1(isAllowed, bool(const V3D &)); + }; + + class MockHKLFilterUnaryLogicOperation : public HKLFilterUnaryLogicOperation { + public: + MockHKLFilterUnaryLogicOperation(const HKLFilter_const_sptr &filter) + : HKLFilterUnaryLogicOperation(filter) {} + + MOCK_CONST_METHOD0(getDescription, std::string()); + MOCK_CONST_METHOD1(isAllowed, bool(const V3D &)); + }; + + class MockHKLFilterBinaryLogicOperation + : public HKLFilterBinaryLogicOperation { + public: + MockHKLFilterBinaryLogicOperation(const HKLFilter_const_sptr &lhs, + const HKLFilter_const_sptr &rhs) + : HKLFilterBinaryLogicOperation(lhs, rhs) {} + + MOCK_CONST_METHOD0(getDescription, std::string()); + MOCK_CONST_METHOD1(isAllowed, bool(const V3D &)); + }; +}; +#endif /* MANTID_GEOMETRY_HKLFILTERTEST_H_ */ diff --git a/Framework/Geometry/test/HKLGeneratorTest.h b/Framework/Geometry/test/HKLGeneratorTest.h new file mode 100644 index 0000000000000000000000000000000000000000..f6cffd236d46002b8dac782b7defd1f04d12bc89 --- /dev/null +++ b/Framework/Geometry/test/HKLGeneratorTest.h @@ -0,0 +1,98 @@ +#ifndef MANTID_GEOMETRY_HKLGENERATORTEST_H_ +#define MANTID_GEOMETRY_HKLGENERATORTEST_H_ + +#include <cxxtest/TestSuite.h> +#include "MantidGeometry/Crystal/HKLGenerator.h" + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; + +class HKLGeneratorTest : public CxxTest::TestSuite { +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static HKLGeneratorTest *createSuite() { return new HKLGeneratorTest(); } + static void destroySuite(HKLGeneratorTest *suite) { delete suite; } + + void test_HKLGeneratorReturnsCorrectSizeSymmetricInt() { + HKLGenerator gen(2, 2, 2); + + TS_ASSERT_EQUALS(gen.size(), 125); + } + + void test_HKLGeneratorReturnsCorrectSizeSymmetricV3D() { + HKLGenerator gen(V3D(2, 2, 2)); + + TS_ASSERT_EQUALS(gen.size(), 125); + } + + void test_HKLGeneratorReturnsCorrectSizeAsymmetricV3D() { + HKLGenerator gen(V3D(-2, -1, -5), V3D(3, 4, -2)); + + TS_ASSERT_EQUALS(gen.size(), 144); + } + + void test_HKLGeneratorReturnsCorrectSizeOne() { + HKLGenerator gen(V3D(-2, -1, -5), V3D(-2, -1, -5)); + TS_ASSERT_EQUALS(gen.size(), 1); + } + + void test_beginIterator() { + HKLGenerator gen(V3D(-2, -2, -3), V3D(1, 1, 0)); + + HKLGenerator::const_iterator it = gen.begin(); + TS_ASSERT_EQUALS(*it, V3D(-2, -2, -3)); + } + + void test_endIterator() { + HKLGenerator gen(V3D(-2, -2, -3), V3D(1, 1, 0)); + + HKLGenerator::const_iterator it = gen.end(); + TS_ASSERT_EQUALS(*it, V3D(2, -2, -3)); + } + + void test_iterator_dereference() { + HKLGenerator::const_iterator it; + TS_ASSERT_EQUALS(*it, V3D(0, 0, 0)); + } + + void test_comparison() { + HKLGenerator::const_iterator it1(V3D(-2, -3, 1)); + HKLGenerator::const_iterator it2(V3D(-2, -3, 1)); + HKLGenerator::const_iterator it3(V3D(-2, -3, 0)); + + TS_ASSERT_EQUALS(it1, it2); + TS_ASSERT_DIFFERS(it1, it3); + } + + void test_iterator_increment() { + HKLGenerator::const_iterator defaultIterator; + TS_ASSERT_THROWS_NOTHING(++defaultIterator); + } + + void test_iterator_dereference_range() { + HKLGenerator::const_iterator it(V3D(-1, -1, -1), V3D(1, 1, 1)); + TS_ASSERT_EQUALS(*it, V3D(-1, -1, -1)); + ++it; + TS_ASSERT_EQUALS(*it, V3D(-1, -1, 0)); + } + + void test_hkl_range() { + HKLGenerator::const_iterator it(V3D(-1, -1, -1), V3D(1, 1, 1)); + HKLGenerator::const_iterator end(V3D(2, -1, -1)); + + std::vector<V3D> hkls; + std::copy(it, end, std::back_inserter(hkls)); + + TS_ASSERT_EQUALS(hkls.size(), 27); + TS_ASSERT_EQUALS(hkls[0], V3D(-1, -1, -1)); + TS_ASSERT_EQUALS(hkls[1], V3D(-1, -1, 0)); + TS_ASSERT_EQUALS(hkls[2], V3D(-1, -1, 1)); + TS_ASSERT_EQUALS(hkls[3], V3D(-1, 0, -1)); + TS_ASSERT_EQUALS(hkls[26], V3D(1, 1, 1)); + + TS_ASSERT_EQUALS(std::distance(it, end), 27); + } +}; + +#endif /* MANTID_GEOMETRY_HKLGENERATORTEST_H_ */ diff --git a/Framework/Geometry/test/ReflectionGeneratorTest.h b/Framework/Geometry/test/ReflectionGeneratorTest.h new file mode 100644 index 0000000000000000000000000000000000000000..9c3da2883aa0664fbd38315122247b895a5123cf --- /dev/null +++ b/Framework/Geometry/test/ReflectionGeneratorTest.h @@ -0,0 +1,126 @@ +#ifndef MANTID_GEOMETRY_REFLECTIONGENERATORTEST_H_ +#define MANTID_GEOMETRY_REFLECTIONGENERATORTEST_H_ + +#include <cxxtest/TestSuite.h> + +#include "MantidGeometry/Crystal/ReflectionGenerator.h" +#include "MantidGeometry/Crystal/BasicHKLFilters.h" + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; + +class ReflectionGeneratorTest : public CxxTest::TestSuite { +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static ReflectionGeneratorTest *createSuite() { + return new ReflectionGeneratorTest(); + } + static void destroySuite(ReflectionGeneratorTest *suite) { delete suite; } + + void test_getUniqueHKLs() { + double dMin = 0.55; + double dMax = 4.0; + + ReflectionGenerator generator( + CrystalStructure("4.126 4.126 4.126", "P m -3 m", "Si 0 0 0 1.0 0.01"), + ReflectionConditionFilter::Centering); + + TS_ASSERT_THROWS_NOTHING(generator.getUniqueHKLs(dMin, dMax)); + + std::vector<V3D> peaks = generator.getUniqueHKLs(dMin, dMax); + + TS_ASSERT_EQUALS(peaks.size(), 68); + TS_ASSERT_EQUALS(peaks[0], V3D(1, 1, 0)); + TS_ASSERT_EQUALS(peaks[11], V3D(3, 2, 0)); + TS_ASSERT_EQUALS(peaks[67], V3D(7, 2, 1)); + + // make d-value list and check that peaks are within limits + std::vector<double> peaksD = generator.getDValues(peaks); + + std::sort(peaksD.begin(), peaksD.end()); + + TS_ASSERT_LESS_THAN_EQUALS(dMin, peaksD.front()); + TS_ASSERT_LESS_THAN_EQUALS(peaksD.back(), dMax); + } + + void test_getHKLs() { + double dMin = 0.55; + double dMax = 4.0; + + // make a structure with P-1 + ReflectionGenerator generator( + CrystalStructure("4.126 4.126 4.126", "P -1", "Si 0 0 0 1.0 0.01"), + ReflectionConditionFilter::Centering); + + std::vector<V3D> unique = generator.getUniqueHKLs(dMin, dMax); + std::vector<V3D> peaks = generator.getHKLs(dMin, dMax); + + // Because of symmetry -1, each reflection has multiplicity 2. + TS_ASSERT_EQUALS(peaks.size(), 2 * unique.size()); + } + + void test_getDValues() { + std::vector<V3D> hkls; + hkls.push_back(V3D(1, 0, 0)); + hkls.push_back(V3D(0, 1, 0)); + hkls.push_back(V3D(0, 0, 1)); + + ReflectionGenerator generator( + CrystalStructure("2 3 5", "P -1", "Si 0 0 0 1.0 0.01")); + std::vector<double> dValues = generator.getDValues(hkls); + + TS_ASSERT_EQUALS(dValues.size(), hkls.size()); + TS_ASSERT_EQUALS(dValues[0], 2.0); + TS_ASSERT_EQUALS(dValues[1], 3.0); + TS_ASSERT_EQUALS(dValues[2], 5.0); + } + + void test_getUniqueHKLsStructureFactor() { + CrystalStructure si("5.43 5.43 5.43", "F m -3 m", + "Si 0.3 0.3 0.3 1.0 0.05"); + + ReflectionGenerator generator(si, + ReflectionConditionFilter::StructureFactor); + + std::vector<V3D> hklsCentering = generator.getUniqueHKLs( + 0.6, 10.0, boost::make_shared<HKLFilterCentering>(si.centering())); + + std::vector<V3D> hklsStructureFactors = generator.getUniqueHKLs(0.6, 10.0); + + TS_ASSERT_EQUALS(hklsCentering.size(), hklsStructureFactors.size()); + + for (size_t i = 0; i < hklsCentering.size(); ++i) { + TS_ASSERT_EQUALS(hklsCentering[i], hklsStructureFactors[i]); + } + } + + void test_getUniqueHKLsHexagonal() { + ReflectionGenerator generator( + CrystalStructure("3.2094 3.2094 5.2108 90.0 90.0 120.0", "P 63/m m c", + "Mg 1/3 2/3 1/4 1.0 0.005"), + ReflectionConditionFilter::StructureFactor); + + std::vector<V3D> hkls = generator.getUniqueHKLs(0.5, 10.0); + + TS_ASSERT_EQUALS(hkls.size(), 88); + + std::vector<double> dValues = generator.getDValues(hkls); + for (size_t i = 0; i < hkls.size(); ++i) { + TS_ASSERT_LESS_THAN(0.5, dValues[i]); + } + } + + void test_getUniqueHKLsTrigonal() { + ReflectionGenerator generator( + CrystalStructure("4.759355 4.759355 12.99231 90.0 90.0 120.0", "R -3 c", + "Al 0 0 0.35217 1.0 0.005; O 0.69365 0 1/4 1.0 0.005"), + ReflectionConditionFilter::StructureFactor); + + std::vector<V3D> hkls = generator.getUniqueHKLs(0.885, 10.0); + + TS_ASSERT_EQUALS(hkls.size(), 44); + } +}; + +#endif /* MANTID_GEOMETRY_REFLECTIONGENERATORTEST_H_ */ diff --git a/Framework/Geometry/test/StructureFactorCalculatorSummationTest.h b/Framework/Geometry/test/StructureFactorCalculatorSummationTest.h new file mode 100644 index 0000000000000000000000000000000000000000..a1ce3c9670a57b3edd0cd2dd6302d8c27c3a33bf --- /dev/null +++ b/Framework/Geometry/test/StructureFactorCalculatorSummationTest.h @@ -0,0 +1,75 @@ +#ifndef MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORSUMMATIONTEST_H_ +#define MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORSUMMATIONTEST_H_ + +#include <cxxtest/TestSuite.h> + +#include "MantidGeometry/Crystal/StructureFactorCalculatorSummation.h" + +#include "MantidGeometry/Crystal/BraggScattererFactory.h" +#include "MantidGeometry/Crystal/SpaceGroupFactory.h" + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; + +class StructureFactorCalculatorSummationTest : public CxxTest::TestSuite { +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static StructureFactorCalculatorSummationTest *createSuite() { + return new StructureFactorCalculatorSummationTest(); + } + static void destroySuite(StructureFactorCalculatorSummationTest *suite) { + delete suite; + } + + void testEquivalentPositionsAreUsed() { + // Approximate crystal structure of Silicon + CrystalStructure si = getCrystalStructure(); + + StructureFactorCalculatorSummation calculator; + calculator.setCrystalStructure(si); + + // {1 0 0} reflections are not allowed because of F centering + TS_ASSERT_LESS_THAN(calculator.getFSquared(V3D(1, 0, 0)), 1e-9); + + // {2 2 2} is forbidden because of Si on special position + TS_ASSERT_LESS_THAN(calculator.getFSquared(V3D(2, 2, 2)), 1e-9); + + // With space group P-1, those are allowed. + si.setSpaceGroup(SpaceGroupFactory::Instance().createSpaceGroup("P -1")); + + calculator.setCrystalStructure(si); + // {1 0 0} reflections are not allowed because of F centering + TS_ASSERT_LESS_THAN(1e-9, calculator.getFSquared(V3D(1, 0, 0))); + + // {2 2 2} is forbidden because of Si on special position + TS_ASSERT_LESS_THAN(1e-9, calculator.getFSquared(V3D(2, 2, 2))); + } + + void testCreateWithFactory() { + CrystalStructure si = getCrystalStructure(); + + StructureFactorCalculator_sptr calculator = + StructureFactorCalculatorFactory::create< + StructureFactorCalculatorSummation>(si); + + // same reflections as in testEquivalentPositionsAreUsed + TS_ASSERT_LESS_THAN(calculator->getFSquared(V3D(1, 0, 0)), 1e-9); + TS_ASSERT_LESS_THAN(calculator->getFSquared(V3D(2, 2, 2)), 1e-9); + } + +private: + CrystalStructure getCrystalStructure() { + CompositeBraggScatterer_sptr scatterers = CompositeBraggScatterer::create(); + scatterers->addScatterer(BraggScattererFactory::Instance().createScatterer( + "IsotropicAtomBraggScatterer", "Element=Si;Position=[0,0,0];U=0.05")); + + CrystalStructure si( + UnitCell(5.43, 5.43, 5.43), + SpaceGroupFactory::Instance().createSpaceGroup("F d -3 m"), scatterers); + + return si; + } +}; + +#endif /* MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORSUMMATIONTEST_H_ */ diff --git a/Framework/Geometry/test/StructureFactorCalculatorTest.h b/Framework/Geometry/test/StructureFactorCalculatorTest.h new file mode 100644 index 0000000000000000000000000000000000000000..72421449f7132c3b281a0f7873be91431a052d59 --- /dev/null +++ b/Framework/Geometry/test/StructureFactorCalculatorTest.h @@ -0,0 +1,107 @@ +#ifndef MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORTEST_H_ +#define MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORTEST_H_ + +#include <cxxtest/TestSuite.h> +#include <gtest/gtest.h> +#include <gmock/gmock.h> + +#include "MantidGeometry/Crystal/StructureFactorCalculator.h" +#include "MantidGeometry/Crystal/SpaceGroupFactory.h" + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; +using ::testing::Mock; +using ::testing::_; +using ::testing::Return; + +class StructureFactorCalculatorTest : public CxxTest::TestSuite { +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static StructureFactorCalculatorTest *createSuite() { + return new StructureFactorCalculatorTest(); + } + static void destroySuite(StructureFactorCalculatorTest *suite) { + delete suite; + } + + void testCrystalStructureSetHookIsCalled() { + CrystalStructure cs(UnitCell(1, 2, 3), + SpaceGroupFactory::Instance().createSpaceGroup("P -1"), + CompositeBraggScatterer::create()); + + MockStructureFactorCalculator calculator; + EXPECT_CALL(calculator, crystalStructureSetHook(_)).Times(1); + + calculator.setCrystalStructure(cs); + + TS_ASSERT(Mock::VerifyAndClearExpectations(&calculator)); + } + + void testGetFSquared() { + MockStructureFactorCalculator calculator; + + EXPECT_CALL(calculator, getF(_)) + .WillRepeatedly(Return(StructureFactor(2.21, 3.1))); + + // Check that the square of 2.21 + i * 3.1 is returned + TS_ASSERT_DELTA(calculator.getFSquared(V3D()), 14.4941, 1e-15); + + TS_ASSERT(Mock::VerifyAndClearExpectations(&calculator)); + } + + void testGetFs() { + MockStructureFactorCalculator calculator; + + int numHKLs = 10; + EXPECT_CALL(calculator, getF(_)) + .Times(numHKLs) + .WillRepeatedly(Return(StructureFactor(2.0, 2.0))); + + std::vector<V3D> hkls(numHKLs); + std::vector<StructureFactor> sfs = calculator.getFs(hkls); + + TS_ASSERT_EQUALS(sfs.size(), hkls.size()); + + for (auto sf = sfs.begin(); sf != sfs.end(); ++sf) { + TS_ASSERT_EQUALS(*sf, StructureFactor(2.0, 2.0)); + } + + TS_ASSERT(Mock::VerifyAndClearExpectations(&calculator)); + } + + void testGetFsSquared() { + MockStructureFactorCalculator calculator; + + int numHKLs = 10; + EXPECT_CALL(calculator, getF(_)) + .Times(numHKLs) + .WillRepeatedly(Return(StructureFactor(2.0, 2.0))); + + std::vector<V3D> hkls(numHKLs); + std::vector<double> sfsSquared = calculator.getFsSquared(hkls); + + TS_ASSERT_EQUALS(sfsSquared.size(), hkls.size()); + + for (auto sf = sfsSquared.begin(); sf != sfsSquared.end(); ++sf) { + TS_ASSERT_EQUALS(*sf, 8.0); + } + + TS_ASSERT(Mock::VerifyAndClearExpectations(&calculator)); + } + +private: + /* This MockStructureFactorCalculator helps to test whether the + * default implementations of the virtual methods work correctly. + * Furthermore it is used to confirm that crystalStructureSetHook + * is called appropriately. + */ + class MockStructureFactorCalculator : public StructureFactorCalculator { + public: + MOCK_CONST_METHOD1(getF, StructureFactor(const V3D &hkl)); + MOCK_METHOD1(crystalStructureSetHook, + void(const CrystalStructure &crystalStructure)); + }; +}; + +#endif /* MANTID_GEOMETRY_STRUCTUREFACTORCALCULATORTEST_H_ */ diff --git a/Framework/PythonInterface/mantid/geometry/src/Exports/CrystalStructure.cpp b/Framework/PythonInterface/mantid/geometry/src/Exports/CrystalStructure.cpp new file mode 100644 index 0000000000000000000000000000000000000000..daf4587a81da370e1ffd0adc07d04f7b625a99c4 --- /dev/null +++ b/Framework/PythonInterface/mantid/geometry/src/Exports/CrystalStructure.cpp @@ -0,0 +1,38 @@ +#include "MantidGeometry/Crystal/CrystalStructure.h" +#include "MantidGeometry/Crystal/IsotropicAtomBraggScatterer.h" +#include <boost/python/class.hpp> +#include <boost/python/make_constructor.hpp> +#include <boost/python/register_ptr_to_python.hpp> + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; +using namespace boost::python; + +namespace { +SpaceGroup_sptr getSpaceGroup(CrystalStructure &self) { + return boost::const_pointer_cast<SpaceGroup>(self.spaceGroup()); +} + +std::vector<std::string> getScatterers(CrystalStructure &self) { + CompositeBraggScatterer_sptr scatterers = self.getScatterers(); + + std::vector<std::string> scattererStrings; + scattererStrings.reserve(scatterers->nScatterers()); + + for (size_t i = 0; i < scatterers->nScatterers(); ++i) { + scattererStrings.push_back( + getIsotropicAtomBraggScattererString(scatterers->getScatterer(i))); + } + + return scattererStrings; +} +} + +void export_CrystalStructure() { + class_<CrystalStructure>("CrystalStructure", no_init) + .def(init<const std::string &, const std::string &, const std::string &>( + (arg("unitCell"), arg("spaceGroup"), arg("scatterers")))) + .def("getUnitCell", &CrystalStructure::cell) + .def("getSpaceGroup", &getSpaceGroup) + .def("getScatterers", &getScatterers); +} diff --git a/Framework/PythonInterface/mantid/geometry/src/Exports/ReflectionGenerator.cpp b/Framework/PythonInterface/mantid/geometry/src/Exports/ReflectionGenerator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..83882c117af6c3965b0e7ed3aae8a75d9dab063b --- /dev/null +++ b/Framework/PythonInterface/mantid/geometry/src/Exports/ReflectionGenerator.cpp @@ -0,0 +1,77 @@ +#include "MantidGeometry/Crystal/ReflectionGenerator.h" +#include "MantidPythonInterface/kernel/Converters/PySequenceToVector.h" +#include <boost/python/class.hpp> +#include <boost/python/enum.hpp> +#include <boost/python/list.hpp> + +using namespace Mantid::Geometry; +using namespace Mantid::Kernel; +using namespace Mantid::PythonInterface; +using namespace boost::python; + +namespace { +boost::python::list getListFromV3DVector(const std::vector<V3D> &hkls) { + boost::python::list hklList; + for (auto it = hkls.begin(); it != hkls.end(); ++it) { + hklList.append(*it); + } + return hklList; +} + +boost::python::list getHKLsDefaultFilter(ReflectionGenerator &self, double dMin, + double dMax) { + return getListFromV3DVector(self.getHKLs(dMin, dMax)); +} + +boost::python::list getHKLsUsingFilter(ReflectionGenerator &self, double dMin, + double dMax, + ReflectionConditionFilter filter) { + return getListFromV3DVector( + self.getHKLs(dMin, dMax, self.getReflectionConditionFilter(filter))); +} + +boost::python::list getUniqueHKLsDefaultFilter(ReflectionGenerator &self, + double dMin, double dMax) { + return getListFromV3DVector(self.getUniqueHKLs(dMin, dMax)); +} + +boost::python::list getUniqueHKLsUsingFilter(ReflectionGenerator &self, + double dMin, double dMax, + ReflectionConditionFilter filter) { + return getListFromV3DVector(self.getUniqueHKLs( + dMin, dMax, self.getReflectionConditionFilter(filter))); +} + +std::vector<double> getDValues(ReflectionGenerator &self, + const boost::python::object &hkls) { + Converters::PySequenceToVector<V3D> converter(hkls); + + return self.getDValues(converter()); +} + +std::vector<double> getFsSquared(ReflectionGenerator &self, + const boost::python::object &hkls) { + Converters::PySequenceToVector<V3D> converter(hkls); + + return self.getFsSquared(converter()); +} +} + +void export_ReflectionGenerator() { + enum_<ReflectionConditionFilter>("ReflectionConditionFilter") + .value("None", ReflectionConditionFilter::None) + .value("Centering", ReflectionConditionFilter::Centering) + .value("SpaceGroup", ReflectionConditionFilter::SpaceGroup) + .value("StructureFactor", ReflectionConditionFilter::StructureFactor) + .export_values(); + + class_<ReflectionGenerator>("ReflectionGenerator", no_init) + .def(init<const CrystalStructure &, optional<ReflectionConditionFilter>>( + (arg("crystalStructure"), arg("defaultFilter")))) + .def("getHKLs", &getHKLsDefaultFilter) + .def("getHKLsUsingFilter", &getHKLsUsingFilter) + .def("getUniqueHKLs", &getUniqueHKLsDefaultFilter) + .def("getUniqueHKLsUsingFilter", &getUniqueHKLsUsingFilter) + .def("getDValues", &getDValues) + .def("getFsSquared", &getFsSquared); +} diff --git a/Framework/PythonInterface/test/python/mantid/geometry/CrystalStructureTest.py b/Framework/PythonInterface/test/python/mantid/geometry/CrystalStructureTest.py new file mode 100644 index 0000000000000000000000000000000000000000..147d8e1e232adbea44574f5a908e45e4a68614ba --- /dev/null +++ b/Framework/PythonInterface/test/python/mantid/geometry/CrystalStructureTest.py @@ -0,0 +1,59 @@ +# pylint: disable=no-init,invalid-name,too-many-public-methods,broad-except +import unittest +from mantid.geometry import CrystalStructure + + +class CrystalStructureTest(unittest.TestCase): + def test_creation(self): + # Some valid constructions + self.assertTrue(self.createCrystalStructureOrRaise("5.43 5.43 5.43", "F d -3 m", "Al 1/3 0.454 1/12 1.0 0.01")) + self.assertTrue(self.createCrystalStructureOrRaise("5.43 5.43 5.43", "C m m m", "Al 1/3 0.454 1/12 1.0 0.01;\n" + "Si 2/3 0.121 1/8")) + self.assertTrue( + self.createCrystalStructureOrRaise("5.43 5.43 5.43 90 90 120", "R -3 c", "Al 1/3 0.454 1/12 1.0 0.01;\n" + "Si 2/3 0.121 1/8")) + + # Invalid unit cell specification + self.assertFalse( + self.createCrystalStructureOrRaise("5.43 5.43 5.43 90.0", "C m m m", "Al 1/3 0.454 1/12 1.0 0.01")) + + # Invalid space group + self.assertFalse( + self.createCrystalStructureOrRaise("5.43 5.43 5.43", "INVALID", "Al 1/3 0.454 1/12 1.0 0.01")) + + # Invalid atom specification + self.assertFalse( + self.createCrystalStructureOrRaise("5.43 5.43 5.43", "C m c e", "Al 1/3 0")) + + def createCrystalStructureOrRaise(self, unitCell, spaceGroup, atomStrings): + try: + CrystalStructure(unitCell, spaceGroup, atomStrings) + return True + except Exception: + return False + + def test_UnitCell(self): + structure = CrystalStructure("5.43 5.42 5.41", "F d -3 m", "Al 1/3 0.454 1/12 1.0 0.01") + cell = structure.getUnitCell() + + self.assertEqual(cell.a(), 5.43) + self.assertEqual(cell.b(), 5.42) + self.assertEqual(cell.c(), 5.41) + + def test_SpaceGroup(self): + structure = CrystalStructure("5.43 5.42 5.41", "F d -3 m", "Al 1/3 0.454 1/12 1.0 0.01") + spaceGroup = structure.getSpaceGroup() + + self.assertEqual(spaceGroup.getHMSymbol(), "F d -3 m") + + def test_scatterers(self): + initialString = "Al 1/3 0.454 1/12 1 0.01;Si 0.1 0.2 0.3 0.99 0.1" + + structure = CrystalStructure("5.43 5.42 5.41", "F d -3 m", initialString) + scatterers = structure.getScatterers() + + self.assertEqual(';'.join(scatterers), initialString) + + +if __name__ == '__main__': + unittest.main() diff --git a/Framework/PythonInterface/test/python/mantid/geometry/ReflectionGeneratorTest.py b/Framework/PythonInterface/test/python/mantid/geometry/ReflectionGeneratorTest.py new file mode 100644 index 0000000000000000000000000000000000000000..4c942f8a7ccd0140d98864097dc3f36f1cfc7ff8 --- /dev/null +++ b/Framework/PythonInterface/test/python/mantid/geometry/ReflectionGeneratorTest.py @@ -0,0 +1,69 @@ +# pylint: disable=no-init,invalid-name,too-many-public-methods +import unittest +from mantid.geometry import CrystalStructure, ReflectionGenerator, ReflectionConditionFilter +from mantid.kernel import V3D +import numpy as np + + +class CrystalStructureTest(unittest.TestCase): + crystalStructure = CrystalStructure("5.431 5.431 5.431", "F d -3 m", "Si 0 0 0 1.0 0.02") + + def test_create(self): + generator = ReflectionGenerator(self.crystalStructure) + generator = ReflectionGenerator(self.crystalStructure, ReflectionConditionFilter.Centering) + + def test_getHKLs(self): + generator = ReflectionGenerator(self.crystalStructure) + hkls = generator.getHKLs(1.0, 10.0) + + self.assertEqual(len(hkls), 138) + + # default is filtering by space group reflection condition + self.assertTrue(V3D(2, 2, 2) in hkls) + self.assertFalse(V3D(1, 0, 0) in hkls) + + def test_getHKLsUsingFilter(self): + generator = ReflectionGenerator(self.crystalStructure) + hkls = generator.getHKLsUsingFilter(1.0, 10.0, ReflectionConditionFilter.StructureFactor) + + self.assertEqual(len(hkls), 130) + + # The 222 is gone now + self.assertFalse(V3D(2, 2, 2) in hkls) + + def test_getUniqueHKLs(self): + generator = ReflectionGenerator(self.crystalStructure) + hkls = generator.getUniqueHKLs(1.0, 10.0) + + self.assertEqual(len(hkls), 9) + + self.assertTrue(V3D(2, 2, 2) in hkls) + self.assertFalse(V3D(1, 0, 0) in hkls) + + def test_getUniqueHKLsUsingFilter(self): + generator = ReflectionGenerator(self.crystalStructure) + hkls = generator.getUniqueHKLsUsingFilter(1.0, 10.0, ReflectionConditionFilter.StructureFactor) + + self.assertEqual(len(hkls), 8) + self.assertFalse(V3D(2, 2, 2) in hkls) + + def test_getDValues(self): + generator = ReflectionGenerator(self.crystalStructure) + hkls = [V3D(1, 0, 0), V3D(1, 1, 1)] + + dValues = generator.getDValues(hkls) + + self.assertEqual(len(hkls), len(dValues)) + self.assertAlmostEqual(dValues[0], 5.431, places=10) + self.assertAlmostEqual(dValues[1], 5.431 / np.sqrt(3.), places=10) + + def test_getFsSquared(self): + generator = ReflectionGenerator(self.crystalStructure) + hkls = generator.getUniqueHKLs(1.0, 10.0) + + fsSquared = generator.getFsSquared(hkls) + + self.assertEqual(len(fsSquared), len(hkls)) + +if __name__ == '__main__': + unittest.main() diff --git a/docs/source/concepts/CrystalStructureAndReflections.rst b/docs/source/concepts/CrystalStructureAndReflections.rst new file mode 100644 index 0000000000000000000000000000000000000000..443dd594566083f6cdd9d048bc7a2e6e3d1b698d --- /dev/null +++ b/docs/source/concepts/CrystalStructureAndReflections.rst @@ -0,0 +1,293 @@ +.. _Crystal structure and reflections: + +Crystal structure and reflections +================================= + +This document describes how crystal structure data can be processed and used in Mantid. For the understanding of the +concepts of :ref:`symmetry <Symmetry groups>` and :ref:`space groups <Point and space groups>` in Mantid it may be +useful to read those introductory articles before proceeding with this document. While there is a short introduction +into theoretical aspects this page is not a replacement for proper introductory text books on the subject where all +these principles are explained in great detail and on a much more general basis. + +Theoretical aspects +~~~~~~~~~~~~~~~~~~~ + +Crystal structures +------------------ + +A crystal is modelled as an infinitely repeating three-dimensional arrangement of scatterers, usually atoms. Due to +the periodic nature of crystals, the entire object can be described by specifying the repeated unit and how +it is repeated. These information are called "crystal structure" and comprise three components: + + 1. Lattice (describes the periodicity of the structure) + 2. Basis (distribution of scatterers in the unit cell) + 3. Space group (describes the symmetry of the arrangement of scatterers) + +The description of the basis depends on the model's level of detail. In the simplest case it could be a list of +point scatterers that are fully characterized by three coordinates (x, y and z in terms of the lattice vectors) and +a scattering length. In reality however, the scatterers are usually atoms that fluctuate about their average position +due to thermal agitation. A basic way to model this motion is to assume it to be an isotropic phenomenon, allowing the +description with one single parameter that quantifies the radius of a sphere in which the scatterer oscillates +harmonically. This so called Debye-Waller-factor will be introduced later on. + +Another important parameter for a basic description of the basis is the so called occupancy. It describes the fraction +of the total number of scatterer-positions that is actually occupied. A common example where this is required are +unordered binary metal mixtures where the same position in the crystal structure is partly filled with two different +atoms in a randomly distributed manner. + +To summarize, a very basic model of the basis comprises a list of scatterers that are in turn characterized by +six parameters: + + 1. Scattering length (known for each element) + 2. Fractional coordinate x + 3. Fractional coordinate y + 4. Fractional coordinate z + 5. Occupancy of the site + 6. Isotropic thermal displacement parameter + +Knowledge of the space group makes it possible to derive all scatterers in the entire unit cell. The symmetry operations +of the space group map each scatterer-position onto a set of equivalent positions that are consequently occupied by the +same type of scatterer as well. Since the unit cell is repeated infinitely in all three directions of space, it is +enough to describe one unit cell. Finally, the lattice is described by six parameters as well, the lengths of the three +lattice vectors (usually given in :math:`\mathrm{\AA{}}`) and the three angles (in degree) between these vectors. + +Reflections and structure factors +--------------------------------- + +In a diffraction experiment the periodic arrangement of atoms is probed using radiation, in this case in the form of +neutrons, of an appropriate wavelength (on the same scale of interatomic distances, typically between 0.5 and +5 :math:`\mathrm{\AA{}}`). The incident beam interacts with the scatterers and in certain orientations the beam is +"reflected" by a flock of lattice planes, a phenomenon which is described by Bragg's law: + +.. math:: + 2d\sin\theta = \lambda + +In this equation :math: `d` is the spacing between the lattice planes, :math:`\theta` is the angle of the incoming beam +and the lattice plane normal and lambda is the wavelength of the radiation. In an experiment theta and lambda are +usually limited, thus they are limiting the range of interplanar spacings that can be probed. In Bragg's law the lattice +plane families are only described by one parameter, the interplanar distance. But each lattice plane family also has an +orientation in space which can be described by the plane's normal vector. Usually the vector is given in terms of the +reciprocal lattice of the structure, where it is reduced to three integers H, K, L, the so called Miller indices. With +knowledge of the :ref:`unit cell <Lattice>` (and thus the :math:`\mathbf{B}`-matrix), the interplanar spacing can also +be computed like this: + +.. math:: + d = \frac{1}{\left|\mathbf{B}\cdot\mathbf{h}\right|} + +The parameters taken into account so far determine the geometric characteristics of Bragg-reflections, i.e. their +position on a detector and their time of flight. But besides these, each reflection also has an intensity. The +intensity is proportional to the squared structure factor, which depends on the kind and arrangement of scatterers in +the unit cell. The structure factor is a complex number and can be calculated for a certain HKL by summing the +contributions of all N atoms j in the unit cell: + +.. math:: + F_{\mathbf{h}} = \sum\limits_{j}^{N}b_j\exp\left(2\pi i \mathbf{h} \cdot \mathbf{x}_j\right) + +In the above equation :math:`b` is the scattering length, :math:`\mathbf{h}` is the Miller index triplet HKL and +:math:`\mathbf{x}` contains the fractional coordinates of the j-th atom. To take into account isotropic thermal +motion of atoms, the term is multiplied with the Debye-Waller factor: + +.. math:: + F_{\mathbf{h}} = \sum\limits_{j}^{N}b_j\exp\left(2\pi i \mathbf{h} \cdot \mathbf{x}_j\right) + \exp\left(-2\pi^2 U/d_{\mathbf{h}}^2\right) + +Here, :math:`U` is the isotropic atomic displacement parameter, usually given in :math:`\mathrm{\AA{}}^2` and +:math:`d` is the lattice spacing discussed above. There are other, more complex models to describe the movement of +atoms, taking into account anisotropic movement and also anharmonic effects. + +Implementation in Mantid +~~~~~~~~~~~~~~~~~~~~~~~~ + +The concepts described above are available through the Python interface of Mantid. Crystal structures are represented +by a class that stores the three necessary pieces of information. Objects of that class can be created by supplying +string representations of those three arguments. + +.. testcode:: ExCrystalStructureConstruction + + from mantid.geometry import CrystalStructure + + silicon = CrystalStructure("5.431 5.431 5.431", "F d -3 m", "Si 0 0 0 1.0 0.05") + + unitCell = silicon.getUnitCell() + print 'Crystal structure of silicon:' + print ' Unit cell:', unitCell.a(), unitCell.b(), unitCell.c(), unitCell.alpha(), unitCell.beta(), unitCell.gamma() + + spaceGroup = silicon.getSpaceGroup() + print ' Space group:', spaceGroup.getHMSymbol() + print ' Point group:', spaceGroup.getPointGroup().getHMSymbol() + + scatterers = silicon.getScatterers() + print ' Total number of scatterers:', len(scatterers) + + for i, scatterer in enumerate(scatterers): + print ' ' + str(i) + ':', scatterer + +The above script produces the following output: + +.. testoutput:: ExCrystalStructureConstruction + + Crystal structure of silicon: + Unit cell: 5.431 5.431 5.431 90.0 90.0 90.0 + Space group: F d -3 m + Point group: m-3m + Total number of scatterers: 1 + 0: Si 0 0 0 1 0.05 + +In general, the unit cell must be specified using either 3 or 6 space-separated floating point numbers, representing +the three axis lengths and the three angles between them. The list of scatterers is required to be a semi-colon +separated list of strings which contain the following information: Element symbol, x, y, z (fractional coordinates), +occupancy (between 0 and 1) and isotropic atomic displacement parameter. The fractional coordinates can also be given +as fractions (for example :math:`1/2` or :math:`1/3`) and for giving the coordinates in hexagonal or trigonal structures +this is highly recommended as there may be precision problems with decimal numbers. + +While the CrystalStructure class is storing information, there is another class that makes use of these information to +generate reflections and calculate structure factors. This class is called ReflectionGenerator and can be constructed +from a CrystalStructure-object: + +.. testcode:: ExReflectionGeneratorConstruction + + from mantid.geometry import CrystalStructure, ReflectionGenerator + from mantid.kernel import V3D + + silicon = CrystalStructure("5.431 5.431 5.431", "F d -3 m", "Si 0 0 0 1.0 0.05") + generator = ReflectionGenerator(silicon) + + # Create list of unique reflections between 0.7 and 3.0 Angstrom + hkls = generator.getUniqueHKLs(0.7, 3.0) + + print 'There are', len(hkls), 'unique reflections for Si in the specified resolution range.' + print 'The reflection [222] is' + (' not' if not V3D(2, 2, 2) in hkls else '') + ' contained in the list.' + +.. testoutput:: ExReflectionGeneratorConstruction + + There are 20 unique reflections for Si in the specified resolution range. + The reflection [222] is contained in the list. + +Checking the reflection conditions of space group :math:`Fd\bar{3}m` (origin choice 1) in the International Tables for +Crystallography shows that if an atom is on the 8a position, additional conditions apply (:math:`h=2n+1` or +:math:`h+k+l=4n` for general reflections). Using these additional conditions, the 222 reflection should in fact +not be in the list. This can be verified by calculating structure factors for the list of reflections and check if +there are very small values present. + +.. testcode:: ExReflectionGeneratorViolations + + from mantid.geometry import CrystalStructure, ReflectionGenerator + import numpy as np + + silicon = CrystalStructure("5.431 5.431 5.431", "F d -3 m", "Si 0 0 0 1.0 0.05") + generator = ReflectionGenerator(silicon) + + # Create list of unique reflections between 0.7 and 3.0 Angstrom + hkls = generator.getUniqueHKLs(0.7, 3.0) + + # Calculate structure factors for those HKLs + fSquared = generator.getFsSquared(hkls) + + # Find HKLs with very small structure factors: + zeroFSquared = [(hkl, sf) for hkl, sf in zip(hkls, fSquared) if sf < 1e-9] + + print 'HKL F^2' + for hkl, sf in zeroFSquared: + print hkl, ' ', np.round(sf, 2) + +The output of the above script should show three reflections with very small values for :math:`F^2`. Their indices +violate the special conditions mentioned in the previous paragraph, so the reflections are actually extinct: + +.. testoutput:: ExReflectionGeneratorViolations + + HKL F^2 + [2,2,2] 0.0 + [4,4,2] 0.0 + [6,2,2] 0.0 + +Those three reflections are included in the list of unique HKLs, because the standard method to determine whether a +reflection is allowed or not uses the space group symmetry which only reflects the general conditions listed in ITA. +It is however possible to exclude those reflections at the cost of more computations by making use of the structure +factor calculation. This can either be done by passing an additional enum-value of the type ReflectionConditionFilter +to the constructor of ReflectionGenerator or by passing it to the actual generator function: + +.. testcode:: ExReflectionGeneratorSF + + from mantid.geometry import CrystalStructure, ReflectionGenerator, ReflectionConditionFilter + from mantid.kernel import V3D + + silicon = CrystalStructure("5.431 5.431 5.431", "F d -3 m", "Si 0 0 0 1.0 0.05") + generator = ReflectionGenerator(silicon) + + # Create list of unique reflections between 0.7 and 3.0 Angstrom, use structure factors for filtering + hkls = generator.getUniqueHKLsUsingFilter(0.7, 3.0, ReflectionConditionFilter.StructureFactor) + + print 'There are', len(hkls), 'unique reflections for Si in the specified resolution range.' + print 'The reflection [222] is' + (' not' if not V3D(2, 2, 2) in hkls else '') + ' contained in the list.' + +With this option, the three reflections from the example above are missing and as an indicator, the [222] reflection +is actually checked: + +.. testoutput:: ExReflectionGeneratorSF + + There are 17 unique reflections for Si in the specified resolution range. + The reflection [222] is not contained in the list. + +Other options for filtering are Centering and None. If the latter one is used the reflections are only filtered +according to their :math:`d`-value to fit the specified range. + +Another capability of ReflectionGenerator is the calculation of :math:`d`-values for a list of HKLs, very similar +to the process for :math:`F^2`: + +.. testcode:: ExReflectionGeneratorCalculateD + + from mantid.geometry import CrystalStructure, ReflectionGenerator, ReflectionConditionFilter + import numpy as np + + silicon = CrystalStructure("5.431 5.431 5.431", "F d -3 m", "Si 0 0 0 1.0 0.05") + generator = ReflectionGenerator(silicon) + + # Create list of unique reflections between 0.7 and 3.0 Angstrom + hkls = generator.getUniqueHKLsUsingFilter(0.7, 3.0, ReflectionConditionFilter.StructureFactor) + + # Calculate d and F^2 + dValues = generator.getDValues(hkls) + fSquared = generator.getFsSquared(hkls) + + pg = silicon.getSpaceGroup().getPointGroup() + + # Make list of tuples and sort by d-values, descending, include point group for multiplicity. + reflections = sorted([(hkl, d, fsq, len(pg.getEquivalents(hkl))) for hkl, d, fsq in zip(hkls, dValues, fSquared)], + key=lambda x: x[1], reverse=True) + + print '{0:<8}{1:>8}{2:>8}{3:>4}'.format('HKL', 'd', 'F^2', 'M') + for reflection in reflections: + print '{0!s:<8}{1:>8.5f}{2:>8.2f}{3:>4}'.format(*reflection) + +This script will print a table with the reflections including their :math:`d`-value, :math:`F^2` and multiplicity due to point group +symmetry: + +.. testoutput:: ExReflectionGeneratorCalculateD + + HKL d F^2 M + [2,2,0] 1.92015 645.02 12 + [3,1,1] 1.63751 263.85 24 + [4,0,0] 1.35775 377.63 6 + [3,3,1] 1.24596 154.47 24 + [4,2,2] 1.10860 221.08 24 + [3,3,3] 1.04520 90.43 8 + [5,1,1] 1.04520 90.43 24 + [4,4,0] 0.96007 129.43 12 + [5,3,1] 0.91801 52.94 48 + [6,2,0] 0.85872 75.78 24 + [5,3,3] 0.82822 31.00 24 + [4,4,4] 0.78390 44.36 8 + [5,5,1] 0.76049 18.15 24 + [7,1,1] 0.76049 18.15 24 + [6,4,2] 0.72575 25.97 48 + [5,5,3] 0.70706 10.62 24 + [7,3,1] 0.70706 10.62 48 + +Further reading +~~~~~~~~~~~~~~~ + +This concept page explains what's available in the Python interface. Some underlying parts may be interesting for C++ +developers, as the concepts of generating and filtering HKLs are pretty much hidden behind the ReflectionGenerator class +in the Python interface. More detail is available in the generated C++ documentation. + +.. categories:: Concepts