UnitCell.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/UnitCell.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/V3D.h"
 
#include <algorithm>
#include <cfloat>
#include <iomanip>
#include <ios>
#include <stdexcept>
 
#include <boost/lexical_cast.hpp>
 
namespace Mantid::Geometry {
using Mantid::Kernel::DblMatrix;
using Mantid::Kernel::V3D;
 
UnitCell::UnitCell() : da(6), ra(6), errorda(6), G(3, 3), Gstar(3, 3), B(3, 3), ModHKL(3, 3), errorModHKL(3, 3) {
  da[0] = da[1] = da[2] = 1.;
  da[3] = da[4] = da[5] = M_PI_2;
  errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] = 0.0;
  MaxOrder = 0;
  CrossTerm = false;
  recalculate();
}
 
UnitCell::UnitCell(double _a, double _b, double _c)
    : da(6), ra(6), errorda(6), G(3, 3), Gstar(3, 3), B(3, 3), ModHKL(3, 3), errorModHKL(3, 3) {
  da[0] = _a;
  da[1] = _b;
  da[2] = _c;
  // Angles are 90 degrees in radians ->Pi/2
  da[3] = da[4] = da[5] = M_PI_2;
  errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] = 0.0;
  MaxOrder = 0;
  CrossTerm = false;
  recalculate();
}
 
UnitCell::UnitCell(double _a, double _b, double _c, double _alpha, double _beta, double _gamma, const int angleunit)
    : da(6), ra(6), errorda(6), G(3, 3), Gstar(3, 3), B(3, 3), ModHKL(3, 3), errorModHKL(3, 3) {
  da[0] = _a;
  da[1] = _b;
  da[2] = _c;
  // Angle transformed in radians
  if (angleunit == angDegrees) {
    da[3] = deg2rad * _alpha;
    da[4] = deg2rad * _beta;
    da[5] = deg2rad * _gamma;
  } else {
    da[3] = _alpha;
    da[4] = _beta;
    da[5] = _gamma;
  }
  errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] = 0.0;
  MaxOrder = 0;
  CrossTerm = false;
  recalculate();
}
 
double UnitCell::a1() const { return da[0]; }
 
double UnitCell::a2() const { return da[1]; }
 
double UnitCell::a3() const { return da[2]; }
double UnitCell::a(int nd) const {
  if (nd < 0 || nd > 2)
    throw(std::invalid_argument("lattice parameter index can change from 0 to 2 "));
  return da[nd];
}
 
double UnitCell::alpha1() const { return da[3]; }
 
double UnitCell::alpha2() const { return da[4]; }
 
double UnitCell::alpha3() const { return da[5]; }
 
double UnitCell::a() const { return da[0]; }
 
double UnitCell::b() const { return da[1]; }
 
double UnitCell::c() const { return da[2]; }
 
double UnitCell::alpha() const { return da[3] * rad2deg; }
 
double UnitCell::beta() const { return da[4] * rad2deg; }
 
double UnitCell::gamma() const { return da[5] * rad2deg; }
 
double UnitCell::b1() const { return ra[0]; }
 
double UnitCell::b2() const { return ra[1]; }
 
double UnitCell::b3() const { return ra[2]; }
 
double UnitCell::beta1() const { return ra[3]; }
 
double UnitCell::beta2() const { return ra[4]; }
 
double UnitCell::beta3() const { return ra[5]; }
 
double UnitCell::astar() const { return ra[0]; }
 
double UnitCell::bstar() const { return ra[1]; }
 
double UnitCell::cstar() const { return ra[2]; }
 
double UnitCell::alphastar() const { return ra[3] * rad2deg; }
 
double UnitCell::betastar() const { return ra[4] * rad2deg; }
 
double UnitCell::gammastar() const { return ra[5] * rad2deg; }
 
double UnitCell::errora() const { return errorda[0]; }
 
double UnitCell::errorb() const { return errorda[1]; }
 
double UnitCell::errorc() const { return errorda[2]; }
 
double UnitCell::erroralpha(const int angleunit) const {
  if (angleunit == angDegrees) {
    return errorda[3] * rad2deg;
  } else {
    return errorda[3];
  }
}
 
double UnitCell::errorbeta(const int angleunit) const {
  if (angleunit == angDegrees) {
    return errorda[4] * rad2deg;
  } else {
    return errorda[4];
  }
}
 
double UnitCell::errorgamma(const int angleunit) const {
  if (angleunit == angDegrees) {
    return errorda[5] * rad2deg;
  } else {
    return errorda[5];
  }
}
 
double UnitCell::errorvolume() const {
  // From latcon.py by Art Schultz
  double V = volume();
  double delta_V_alphaV = 0.0;
  if (erroralpha() > 0.0) {
    double alpha1 = alpha() - 0.5 * erroralpha();
    double Va1 = UnitCell(a(), b(), c(), alpha1, beta(), gamma()).volume();
    double alpha2 = alpha() + 0.5 * erroralpha();
    double Va2 = UnitCell(a(), b(), c(), alpha2, beta(), gamma()).volume();
    delta_V_alphaV = (Va2 - Va1) / V;
  }
 
  double delta_V_betaV = 0.0;
  if (errorbeta() > 0.0) {
    double beta1 = beta() - 0.5 * errorbeta();
    double Va1 = UnitCell(a(), b(), c(), alpha(), beta1, gamma()).volume();
    double beta2 = beta() + 0.5 * errorbeta();
    double Va2 = UnitCell(a(), b(), c(), alpha(), beta2, gamma()).volume();
    delta_V_betaV = (Va2 - Va1) / V;
  }
 
  double delta_V_gammaV = 0.0;
  if (errorgamma() > 0.0) {
    double gamma1 = gamma() - 0.5 * errorgamma();
    double Va1 = UnitCell(a(), b(), c(), alpha(), beta(), gamma1).volume();
    double gamma2 = gamma() + 0.5 * errorgamma();
    double Va2 = UnitCell(a(), b(), c(), alpha(), beta(), gamma2).volume();
    delta_V_gammaV = (Va2 - Va1) / V;
  }
 
  return V * sqrt(std::pow(errora() / a(), 2) + std::pow(errorb() / b(), 2) + std::pow(errorc() / c(), 2) +
                  std::pow(delta_V_alphaV, 2) + std::pow(delta_V_betaV, 2) + std::pow(delta_V_gammaV, 2));
}
 
void UnitCell::set(double _a, double _b, double _c, double _alpha, double _beta, double _gamma, const int angleunit) {
  da[0] = _a;
  da[1] = _b;
  da[2] = _c;
  if (angleunit == angDegrees) {
    da[3] = deg2rad * _alpha;
    da[4] = deg2rad * _beta;
    da[5] = deg2rad * _gamma;
  } else {
    da[3] = _alpha;
    da[4] = _beta;
    da[5] = _gamma;
  }
  recalculate();
}
 
void UnitCell::setError(double _aerr, double _berr, double _cerr, double _alphaerr, double _betaerr, double _gammaerr,
                        const int angleunit) {
  errorda[0] = _aerr;
  errorda[1] = _berr;
  errorda[2] = _cerr;
  if (angleunit == angDegrees) {
    errorda[3] = deg2rad * _alphaerr;
    errorda[4] = deg2rad * _betaerr;
    errorda[5] = deg2rad * _gammaerr;
  } else {
    errorda[3] = _alphaerr;
    errorda[4] = _betaerr;
    errorda[5] = _gammaerr;
  }
}
 
void UnitCell::setModHKL(double _dh1, double _dk1, double _dl1, double _dh2, double _dk2, double _dl2, double _dh3,
                         double _dk3, double _dl3) {
  ModHKL[0][0] = _dh1;
  ModHKL[1][0] = _dk1;
  ModHKL[2][0] = _dl1;
  ModHKL[0][1] = _dh2;
  ModHKL[1][1] = _dk2;
  ModHKL[2][1] = _dl2;
  ModHKL[0][2] = _dh3;
  ModHKL[1][2] = _dk3;
  ModHKL[2][2] = _dl3;
}
 
void UnitCell::setModHKL(const DblMatrix &newModHKL) { ModHKL = newModHKL; }
 
void UnitCell::setErrorModHKL(const DblMatrix &newErrorModHKL) { errorModHKL = newErrorModHKL; }
 
void UnitCell::setErrorModHKL(double _dh1err, double _dk1err, double _dl1err, double _dh2err, double _dk2err,
                              double _dl2err, double _dh3err, double _dk3err, double _dl3err) {
  errorModHKL[0][0] = _dh1err;
  errorModHKL[1][0] = _dk1err;
  errorModHKL[2][0] = _dl1err;
  errorModHKL[0][1] = _dh2err;
  errorModHKL[1][1] = _dk2err;
  errorModHKL[2][1] = _dl2err;
  errorModHKL[0][2] = _dh3err;
  errorModHKL[1][2] = _dk3err;
  errorModHKL[2][2] = _dl3err;
}
 
void UnitCell::setModVec1(double _dh1, double _dk1, double _dl1) {
  ModHKL[0][0] = _dh1;
  ModHKL[1][0] = _dk1;
  ModHKL[2][0] = _dl1;
}
 
void UnitCell::setModVec2(double _dh2, double _dk2, double _dl2) {
  ModHKL[0][1] = _dh2;
  ModHKL[1][1] = _dk2;
  ModHKL[2][1] = _dl2;
}
 
void UnitCell::setModVec3(double _dh3, double _dk3, double _dl3) {
  ModHKL[0][2] = _dh3;
  ModHKL[1][2] = _dk3;
  ModHKL[2][2] = _dl3;
}
 
void UnitCell::setModVec1(const V3D &newModVec) {
  ModHKL[0][0] = newModVec[0];
  ModHKL[1][0] = newModVec[1];
  ModHKL[2][0] = newModVec[2];
}
 
void UnitCell::setModVec2(const V3D &newModVec) {
  ModHKL[0][1] = newModVec[0];
  ModHKL[1][1] = newModVec[1];
  ModHKL[2][1] = newModVec[2];
}
 
void UnitCell::setModVec3(const V3D &newModVec) {
  ModHKL[0][2] = newModVec[0];
  ModHKL[1][2] = newModVec[1];
  ModHKL[2][2] = newModVec[2];
}
 
void UnitCell::setModerr(int i, double _dherr, double _dkerr, double _dlerr) {
  errorModHKL[0][i] = _dherr;
  errorModHKL[1][i] = _dkerr;
  errorModHKL[2][i] = _dlerr;
}
 
void UnitCell::setModerr1(double _dh1err, double _dk1err, double _dl1err) {
  errorModHKL[0][0] = _dh1err;
  errorModHKL[1][0] = _dk1err;
  errorModHKL[2][0] = _dl1err;
}
 
void UnitCell::setModerr2(double _dh2err, double _dk2err, double _dl2err) {
  errorModHKL[0][1] = _dh2err;
  errorModHKL[1][1] = _dk2err;
  errorModHKL[2][1] = _dl2err;
}
 
void UnitCell::setModerr3(double _dh3err, double _dk3err, double _dl3err) {
  errorModHKL[0][2] = _dh3err;
  errorModHKL[1][2] = _dk3err;
  errorModHKL[2][2] = _dl3err;
}
 
void UnitCell::setMaxOrder(int MaxO) { MaxOrder = MaxO; }
 
void UnitCell::setCrossTerm(bool CT) { CrossTerm = CT; }
 
const Kernel::V3D UnitCell::getModVec(int j) const { return V3D(getdh(j), getdk(j), getdl(j)); }
 
const Kernel::V3D UnitCell::getVecErr(int j) const { return V3D(getdherr(j), getdkerr(j), getdlerr(j)); }
 
const Kernel::DblMatrix &UnitCell::getModHKL() const { return ModHKL; }
 
double UnitCell::getdh(int j) const { return ModHKL[0][j]; }
 
double UnitCell::getdk(int j) const { return ModHKL[1][j]; }
 
double UnitCell::getdl(int j) const { return ModHKL[2][j]; }
 
double UnitCell::getdherr(int j) const { return errorModHKL[0][j]; }
 
double UnitCell::getdkerr(int j) const { return errorModHKL[1][j]; }
 
double UnitCell::getdlerr(int j) const { return errorModHKL[2][j]; }
 
int UnitCell::getMaxOrder() const { return MaxOrder; }
 
bool UnitCell::getCrossTerm() const { return CrossTerm; }
 
void UnitCell::seta(double _a) {
  da[0] = _a;
  recalculate();
}
void UnitCell::setErrora(double _aerr) { errorda[0] = _aerr; }
 
void UnitCell::setb(double _b) {
  da[1] = _b;
  recalculate();
}
void UnitCell::setErrorb(double _berr) { errorda[1] = _berr; }
void UnitCell::setc(double _c) {
  da[2] = _c;
  recalculate();
}
void UnitCell::setErrorc(double _cerr) { errorda[2] = _cerr; }
void UnitCell::setalpha(double _alpha, const int angleunit) {
  if (angleunit == angDegrees)
    da[3] = deg2rad * _alpha;
  else
    da[3] = _alpha;
  recalculate();
}
void UnitCell::setErroralpha(double _alphaerr, const int angleunit) {
  if (angleunit == angDegrees)
    errorda[3] = deg2rad * _alphaerr;
  else
    errorda[3] = _alphaerr;
}
void UnitCell::setbeta(double _beta, const int angleunit) {
  if (angleunit == angDegrees)
    da[4] = deg2rad * _beta;
  else
    da[4] = _beta;
  recalculate();
}
 
void UnitCell::setErrorbeta(double _betaerr, const int angleunit) {
  if (angleunit == angDegrees)
    errorda[4] = deg2rad * _betaerr;
  else
    errorda[4] = _betaerr;
}
 
void UnitCell::setgamma(double _gamma, const int angleunit) {
  if (angleunit == angDegrees)
    da[5] = deg2rad * _gamma;
  else
    da[5] = _gamma;
  recalculate();
}
 
void UnitCell::setErrorgamma(double _gammaerr, const int angleunit) {
  if (angleunit == angDegrees)
    errorda[5] = deg2rad * _gammaerr;
  else
    errorda[5] = _gammaerr;
}
 
double UnitCell::d(double h, double k, double l) const { return 1.0 / dstar(V3D(h, k, l)); }
 
double UnitCell::d(const V3D &hkl) const { return 1.0 / dstar(hkl); }
 
double UnitCell::dstar(double h, double k, double l) const {
  return dstar(V3D(h, k, l)); // create a V3D vector h,k,l
}
 
double UnitCell::dstar(const V3D &hkl) const {
  V3D Q = B * hkl; // transform into $AA^-1$
  return Q.norm();
}
 
double UnitCell::recAngle(double h1, double k1, double l1, double h2, double k2, double l2, const int angleunit) const {
  V3D Q1(h1, k1, l1), Q2(h2, k2, l2);
  double E, ang;
  Q1 = Gstar * Q1;
  E = Q1.scalar_prod(Q2);
  double temp = E / dstar(h1, k1, l1) / dstar(h2, k2, l2);
  if (temp > 1)
    ang = 0.;
  else if (temp < -1)
    ang = M_PI;
  else
    ang = acos(temp);
  if (angleunit == angDegrees)
    return rad2deg * ang;
  else
    return ang;
}
 
double UnitCell::volume() const {
  double volume = G.determinant();
  return sqrt(volume);
}
 
double UnitCell::recVolume() const {
  double recvolume = Gstar.determinant();
  return sqrt(recvolume);
}
 
const Kernel::DblMatrix &UnitCell::getG() const { return G; }
 
const Kernel::DblMatrix &UnitCell::getGstar() const { return Gstar; }
 
const Kernel::DblMatrix &UnitCell::getB() const { return B; }
 
const Kernel::DblMatrix &UnitCell::getBinv() const { return Binv; }
 
void UnitCell::recalculate() {
  if ((da[3] > da[4] + da[5]) || (da[4] > da[3] + da[5]) || (da[5] > da[4] + da[3])) {
    throw std::invalid_argument("Invalid angles");
  }
  calculateG();
  calculateGstar();
  calculateReciprocalLattice();
  calculateB();
}
 
void UnitCell::calculateG() {
  G[0][0] = da[0] * da[0];
  G[1][1] = da[1] * da[1];
  G[2][2] = da[2] * da[2];
  G[0][1] = da[0] * da[1] * cos(da[5]);
  G[0][2] = da[0] * da[2] * cos(da[4]);
  G[1][2] = da[1] * da[2] * cos(da[3]);
  G[1][0] = G[0][1];
  G[2][0] = G[0][2];
  G[2][1] = G[1][2];
}
 
void UnitCell::calculateGstar() {
  // Reciprocal metrix tensor is simply the inverse of the direct one
  double det = G.determinant();
  if (det == 0) {
    throw std::range_error("UnitCell not properly initialized");
  }
  Gstar = G;
  if (Gstar.Invert() == 0) {
    throw std::range_error("UnitCell not properly initialized");
  }
}
 
void UnitCell::calculateReciprocalLattice() {
  ra[0] = sqrt(Gstar[0][0]); // a*
  ra[1] = sqrt(Gstar[1][1]); // b*
  ra[2] = sqrt(Gstar[2][2]); // c*
  auto acosThreshold = [](double x) { return std::abs(x) > 1e-15 ? acos(x) : M_PI_2; };
  ra[3] = acosThreshold(Gstar[1][2] / ra[1] / ra[2]); // alpha*
  ra[4] = acosThreshold(Gstar[0][2] / ra[0] / ra[2]); // beta*
  ra[5] = acosThreshold(Gstar[0][1] / ra[0] / ra[1]); // gamma*
}
 
void UnitCell::calculateB() {
  // B matrix using a right handed coordinate system with b1 along x and y in
  // the (b1,b2) plane.
  // This is the convention in Busing and Levy.
  // | b1 b2cos(beta3)      b3cos(beta2)        |
  // | 0  b2sin(beta3) -b3sin(beta2)cos(alpha1) |
  // | 0       0                  1/a3          |
  B[0][0] = ra[0];
  B[0][1] = ra[1] * cos(ra[5]);
  B[0][2] = ra[2] * cos(ra[4]);
  B[1][0] = 0.;
  B[1][1] = ra[1] * sin(ra[5]);
  B[1][2] = -ra[2] * sin(ra[4]) * cos(da[3]);
  B[2][0] = 0.;
  B[2][1] = 0.;
  B[2][2] = 1. / da[2];
 
  Binv = B;
  Binv.Invert();
}
 
void UnitCell::recalculateFromGstar(const DblMatrix &NewGstar) {
  if (NewGstar.numRows() != 3 || NewGstar.numCols() != 3) {
    std::ostringstream msg;
    msg << "UnitCell::recalculateFromGstar - Expected a 3x3 matrix but was "
           "given a "
        << NewGstar.numRows() << "x" << NewGstar.numCols();
    throw std::invalid_argument(msg.str());
  }
 
  if (NewGstar[0][0] * NewGstar[1][1] * NewGstar[2][2] <= 0.)
    throw std::invalid_argument("NewGstar");
  Gstar = NewGstar;
  calculateReciprocalLattice();
  G = Gstar;
  G.Invert();
  da[0] = sqrt(G[0][0]);                 // a
  da[1] = sqrt(G[1][1]);                 // b
  da[2] = sqrt(G[2][2]);                 // c
  da[3] = acos(G[1][2] / da[1] / da[2]); // alpha
  da[4] = acos(G[0][2] / da[0] / da[2]); // beta
  da[5] = acos(G[0][1] / da[0] / da[1]); // gamma
  calculateB();
}
 
bool UnitCell::operator==(const UnitCell &other) const {
  return da == other.da; // da error not used in comparison
}
bool UnitCell::operator!=(const UnitCell &other) const { return !this->operator==(other); }
 
std::ostream &operator<<(std::ostream &out, const UnitCell &unitCell) {
  // always show the lattice constants
  out << "Lattice Parameters:" << std::fixed << std::setprecision(6) << std::setw(12) << unitCell.a() << std::fixed
      << std::setprecision(6) << std::setw(12) << unitCell.b() << std::fixed << std::setprecision(6) << std::setw(12)
      << unitCell.c() << std::fixed << std::setprecision(6) << std::setw(12) << unitCell.alpha() << std::fixed
      << std::setprecision(6) << std::setw(12) << unitCell.beta() << std::fixed << std::setprecision(6) << std::setw(12)
      << unitCell.gamma() << std::fixed << std::setprecision(6) << " " << std::setw(12) << unitCell.volume();
 
  // write out the uncertainty if there is a positive one somewhere
  if ((unitCell.errora() > 0) || (unitCell.errorb() > 0) || (unitCell.errorc() > 0) || (unitCell.erroralpha() > 0) ||
      (unitCell.errorbeta() > 0) || (unitCell.errorgamma() > 0))
    out << "\nParameter Errors  :" << std::fixed << std::setprecision(6) << std::setw(12) << unitCell.errora()
        << std::fixed << std::setprecision(6) << std::setw(12) << unitCell.errorb() << std::fixed
        << std::setprecision(6) << std::setw(12) << unitCell.errorc() << std::fixed << std::setprecision(6)
        << std::setw(12) << unitCell.erroralpha() << std::fixed << std::setprecision(6) << std::setw(12)
        << unitCell.errorbeta() << std::fixed << std::setprecision(6) << std::setw(12) << unitCell.errorgamma()
        << std::fixed << std::setprecision(6) << std::setw(12) << unitCell.errorvolume();
 
  return out;
}
 
std::string unitCellToStr(const UnitCell &unitCell) {
  std::ostringstream stream;
  stream << std::setprecision(9);
 
  stream << unitCell.a() << " " << unitCell.b() << " " << unitCell.c() << " " << unitCell.alpha() << " "
         << unitCell.beta() << " " << unitCell.gamma();
 
  return stream.str();
}
 
UnitCell strToUnitCell(const std::string &unitCellString) {
 
  Mantid::Kernel::StringTokenizer cellTokens(unitCellString, " ", Mantid::Kernel::StringTokenizer::TOK_IGNORE_EMPTY);
 
  std::vector<double> components;
  components.reserve(cellTokens.size());
  std::transform(cellTokens.cbegin(), cellTokens.cend(), std::back_inserter(components),
                 [](const auto &token) { return boost::lexical_cast<double>(token); });
 
  switch (components.size()) {
  case 3:
    return UnitCell(components[0], components[1], components[2]);
  case 6:
    return UnitCell(components[0], components[1], components[2], components[3], components[4], components[5]);
  default:
    throw std::runtime_error("Failed to parse unit cell input string: " + unitCellString);
  }
}
 
} // namespace Mantid::Geometry