SymmetryOperation.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidGeometry/Crystal/SymmetryOperation.h"
#include "MantidGeometry/Crystal/SymmetryOperationSymbolParser.h"
#include <memory>
 
#include <boost/algorithm/string.hpp>
 
namespace Mantid::Geometry {
 
SymmetryOperation::SymmetryOperation()
    : m_order(1), m_transposedInverseMatrix(Kernel::IntMatrix(3, 3, true)), m_reducedVector(), m_identifier(),
      m_matrixVectorPair() {
  m_identifier = SymmetryOperationSymbolParser::getNormalizedIdentifier(m_matrixVectorPair);
}
 
SymmetryOperation::SymmetryOperation(const std::string &identifier) {
  init(SymmetryOperationSymbolParser::parseIdentifier(identifier));
}
 
SymmetryOperation::SymmetryOperation(const Kernel::IntMatrix &matrix, const V3R &vector) {
  init(MatrixVectorPair<int, V3R>(matrix, vector));
}
 
SymmetryOperation::SymmetryOperation(const Mantid::Kernel::DblMatrix &matrix, const Mantid::Geometry::V3R &vector) {
  init(MatrixVectorPair<int, V3R>(convertMatrix<int>(matrix), vector));
}
 
void SymmetryOperation::init(const MatrixVectorPair<int, V3R> &matrixVectorPair) {
  m_matrixVectorPair = matrixVectorPair;
 
  // Inverse matrix for HKL operations.
  m_transposedInverseMatrix = m_matrixVectorPair.getMatrix();
  m_transposedInverseMatrix.Invert();
  m_transposedInverseMatrix = m_transposedInverseMatrix.Transpose();
 
  m_order = getOrderFromMatrix(m_matrixVectorPair.getMatrix());
  m_identifier = SymmetryOperationSymbolParser::getNormalizedIdentifier(m_matrixVectorPair);
 
  m_reducedVector = getReducedVector(m_matrixVectorPair.getMatrix(), m_matrixVectorPair.getVector());
}
 
const Kernel::IntMatrix &SymmetryOperation::matrix() const { return m_matrixVectorPair.getMatrix(); }
 
const V3R &SymmetryOperation::vector() const { return m_matrixVectorPair.getVector(); }
 
const V3R &SymmetryOperation::reducedVector() const { return m_reducedVector; }
 
size_t SymmetryOperation::order() const { return m_order; }
 
std::string SymmetryOperation::identifier() const { return m_identifier; }
 
bool SymmetryOperation::isIdentity() const {
  return !hasTranslation() && m_matrixVectorPair.getMatrix() == Kernel::IntMatrix(3, 3, true);
}
 
bool SymmetryOperation::hasTranslation() const { return m_matrixVectorPair.getVector() != 0; }
 
Kernel::V3D SymmetryOperation::transformHKL(const Kernel::V3D &hkl) const { return m_transposedInverseMatrix * hkl; }
 
SymmetryOperation SymmetryOperation::operator*(const SymmetryOperation &operand) const {
  MatrixVectorPair<int, V3R> result = m_matrixVectorPair * operand.m_matrixVectorPair;
  return SymmetryOperation(result.getMatrix(), result.getVector());
}
 
SymmetryOperation SymmetryOperation::inverse() const {
  MatrixVectorPair<int, V3R> inverse = m_matrixVectorPair.getInverse();
 
  return SymmetryOperation(inverse.getMatrix(), inverse.getVector());
}
 
SymmetryOperation SymmetryOperation::operator^(size_t exponent) const {
  // If the exponent is 1, no calculations are necessary.
  if (exponent == 1) {
    return SymmetryOperation(*this);
  }
 
  SymmetryOperation op;
 
  // The same for 0, which means identity in every case.
  if (exponent == 0) {
    return op;
  }
 
  for (size_t i = 0; i < exponent; ++i) {
    op = (*this) * op;
  }
 
  return op;
}
 
bool SymmetryOperation::operator==(const SymmetryOperation &other) const {
  return m_matrixVectorPair == other.m_matrixVectorPair;
}
 
bool SymmetryOperation::operator<(const SymmetryOperation &other) const { return m_identifier < other.m_identifier; }
 
bool SymmetryOperation::operator!=(const SymmetryOperation &other) const { return !(this->operator==(other)); }
 
size_t SymmetryOperation::getOrderFromMatrix(const Kernel::IntMatrix &matrix) const {
  int trace = matrix.Trace();
  int determinant = matrix.determinant();
 
  if (determinant == 1) {
    switch (trace) {
    case 3:
      return 1;
    case 2:
      return 6;
    case 1:
      return 4;
    case 0:
      return 3;
    case -1:
      return 2;
    default:
      break;
    }
  } else if (determinant == -1) {
    switch (trace) {
    case -3:
      return 2;
    case -2:
      return 6;
    case -1:
      return 4;
    case 0:
      return 6;
    case 1:
      return 2;
    default:
      break;
    }
  }
 
  throw std::runtime_error("There is something wrong with supplied matrix.");
}
 
V3R SymmetryOperation::getReducedVector(const Kernel::IntMatrix &matrix, const V3R &vector) const {
  Kernel::IntMatrix translationMatrix(3, 3, false);
 
  for (size_t i = 0; i < order(); ++i) {
    Kernel::IntMatrix tempMatrix(3, 3, true);
 
    for (size_t j = 0; j < i; ++j) {
      tempMatrix *= matrix;
    }
 
    translationMatrix += tempMatrix;
  }
 
  return (translationMatrix * vector) * RationalNumber(1, static_cast<int>(order()));
}
 
V3R getWrappedVector(const V3R &vector) {
  V3R wrappedVector(vector);
 
  for (size_t i = 0; i < 3; ++i) {
    wrappedVector[i] -= (vector[i].numerator() / vector[i].denominator());
 
    if (wrappedVector[i] < 0) {
      wrappedVector[i] += 1;
    }
  }
 
  return wrappedVector;
}
 
Kernel::V3D getWrappedVector(const Kernel::V3D &vector) {
  Kernel::V3D wrappedVector;
 
  for (size_t i = 0; i < 3; ++i) {
    wrappedVector[i] = fmod(vector[i], 1.0);
 
    if (wrappedVector[i] < 0) {
      wrappedVector[i] += 1.0;
    }
  }
 
  return wrappedVector;
}
 
std::ostream &operator<<(std::ostream &stream, const SymmetryOperation &operation) {
  stream << "[" + operation.identifier() + "]";
 
  return stream;
}
 
std::istream &operator>>(std::istream &stream, SymmetryOperation &operation) {
  std::string identifier;
  std::getline(stream, identifier);
 
  size_t i = identifier.find_first_of('[');
  size_t j = identifier.find_last_of(']');
 
  if (i == std::string::npos || j == std::string::npos) {
    throw std::runtime_error("Cannot construct SymmetryOperation from identifier: " + identifier);
  }
 
  operation = SymmetryOperation(identifier.substr(i + 1, j - i - 1));
 
  return stream;
}
 
SymmetryOperation getUnitCellIntervalOperation(const SymmetryOperation &symOp) {
  return SymmetryOperation(symOp.matrix(), getWrappedVector(symOp.vector()));
}
 
} // namespace Mantid::Geometry