CompositionConverter.cc

Go to the documentation of this file.

#include "casm/clex/CompositionConverter.hh"
 
#include "casm/misc/CASM_Eigen_math.hh"
 
// currently still relies on this for getting standard composition axes
#include "casm/clex/ParamComposition.hh"
#include "casm/misc/algorithm.hh"
 
namespace CASM {
namespace Local {
static bool _is_vacancy(const std::string &name) {
  return (name == "VA" || name == "va" || name == "Va");
}
}  // namespace Local
CompositionConverter::size_type CompositionConverter::independent_compositions()
    const {
  return m_to_x.rows();
}
 
std::string CompositionConverter::comp_var(size_type i) {
  return std::string(1, (char)(i + (int)'a'));
}
 
std::vector<std::string> CompositionConverter::components() const {
  return m_components;
}
 
Eigen::VectorXd CompositionConverter::origin() const { return m_origin; }
 
Eigen::VectorXd CompositionConverter::end_member(size_type i) const {
  return m_end_members.col(i);
}
 
Eigen::MatrixXd CompositionConverter::dparam_dmol() const { return m_to_x; }
 
Eigen::MatrixXd CompositionConverter::dmol_dparam() const { return m_to_n; }
 
Eigen::VectorXd CompositionConverter::param_composition(
    const Eigen::VectorXd &n) const {
  return m_to_x * (n - m_origin);
}
 
Eigen::VectorXd CompositionConverter::dparam_composition(
    const Eigen::VectorXd &dn) const {
  return m_to_x * dn;
}
 
Eigen::VectorXd CompositionConverter::mol_composition(
    const Eigen::VectorXd &x) const {
  return m_origin + m_to_n * x;
}
 
Eigen::VectorXd CompositionConverter::dmol_composition(
    const Eigen::VectorXd &dx) const {
  return m_to_n * dx;
}
 
Eigen::VectorXd CompositionConverter::param_chem_pot(
    const Eigen::VectorXd chem_pot) const {
  return m_to_n.transpose() * chem_pot;
}
 
std::string CompositionConverter::mol_formula() const {
  // n = origin + m_to_n * x
 
  std::stringstream tstr;
 
  // for each molecule:
  for (int i = 0; i < m_components.size(); i++) {
    bool first_char = true;
 
    // print mol name 'A('
    tstr << m_components[i] << "(";
 
    // constant term from origin
    // print 'x' if x != 0
    if (!almost_zero(m_origin(i))) {
      first_char = false;
      tstr << m_origin(i);
    }
 
    // terms from m_to_n columns
    for (int j = 0; j < m_to_n.cols(); j++) {
      // print nothing if x == 0
      if (almost_zero(m_to_n(i, j))) {
        continue;
      }
 
      // print '+' if x>0 mid-expression
      if (!first_char && m_to_n(i, j) > 0) tstr << '+';
      // print '-' if x<0
      else if (m_to_n(i, j) < 0)
        tstr << '-';
      // print absolute value of x if |x|!=1
      if (!almost_equal(std::abs(m_to_n(i, j)), 1.0))
        tstr << std::abs(m_to_n(i, j));
      // print variable ('a','b',etc...)
      tstr << comp_var(j);
 
      first_char = false;
    }
 
    // close ')'
    tstr << ")";
  }
 
  return tstr.str();
}
 
std::string CompositionConverter::param_formula() const {
  // x_i = m_to_x*(n - origin) = m_to_x*n - m_to_x*origin
 
  std::stringstream tstr;
 
  Eigen::VectorXd v = -m_to_x * m_origin;
 
  // for each independent composition:
  for (int i = 0; i < independent_compositions(); i++) {
    bool first_char = true;
 
    // print mol name 'a('
    tstr << comp_var(i) << "(";
 
    // constant term from origin
    // print 'x' if x != 0
    if (!almost_zero(v(i))) {
      first_char = false;
      tstr << v(i);
    }
 
    // terms from m_to_x columns
    for (int j = 0; j < m_to_x.cols(); j++) {
      double coeff = m_to_x(i, j);
 
      // print nothing if n == 0
      if (almost_zero(coeff)) {
        continue;
      }
 
      // print 'A' or '+A' if n == 1
      if (almost_zero(coeff - 1)) {
        if (!first_char) {
          tstr << '+';
        }
        tstr << m_components[j];
      }
 
      // print '-A' if n == -1
      else if (almost_zero(coeff + 1)) {
        tstr << '-' << m_components[j];
      }
 
      // print 'nA' or '+nA' if n > 0
      else if (coeff > 0) {
        if (!first_char) {
          tstr << '+';
        }
        tstr << coeff << m_components[j];
      }
 
      // print '-nA' if n < 0
      else {
        tstr << coeff << m_components[j];
      }
 
      first_char = false;
    }
 
    // close ')'
    tstr << ")";
  }
 
  return tstr.str();
}
 
std::string CompositionConverter::comp_formula(size_type i) const {
  // comp(i) = m_to_x(i,j)*(comp_n(j) - m_origin(j)) + ...
 
  std::stringstream ss;
 
  auto comp_x_str = [&]() { return "comp(" + comp_var(i) + ")"; };
 
  auto comp_n_str = [&](int j) { return "comp_n(" + m_components[j] + ")"; };
 
  auto delta_str = [&](int j) {
    std::stringstream tss;
    // print '(comp_n(J) - m_origin(j))' if m_origin(j) != 0
    if (!almost_zero(m_origin(j))) {
      tss << "(" << comp_n_str(j) << " - " << m_origin(j) << ")";
    }
    // print 'comp_n(J)'
    else {
      tss << comp_n_str(j);
    }
    return tss.str();
  };
 
  ss << comp_x_str() << " = ";
  bool first_term = true;
  for (int j = 0; j < m_to_x.cols(); ++j) {
    double coeff = m_to_x(i, j);
 
    // print nothing if coeff == 0
    if (almost_zero(coeff)) {
      continue;
    }
 
    // if coeff < 0
    if (coeff < 0) {
      if (!first_term) {
        ss << " - " << -coeff << "*" << delta_str(j);
      } else {
        ss << coeff << "*" << delta_str(j);
      }
    }
 
    // if coeff > 0
    else {
      if (!first_term) {
        ss << " + " << coeff << "*" << delta_str(j);
      } else {
        ss << coeff << "*" << delta_str(j);
      }
    }
    ss << " ";
 
    first_term = false;
  }
 
  return ss.str();
}
 
std::string CompositionConverter::comp_n_formula(size_type i) const {
  // comp_n(i) = m_origin(j) + m_to_n(i,j)*comp(j) + ...
 
  std::stringstream ss;
 
  auto comp_x_str = [&](int j) { return "comp(" + comp_var(j) + ")"; };
 
  auto comp_n_str = [&](int j) { return "comp_n(" + m_components[j] + ")"; };
 
  ss << comp_n_str(i) << " = ";
  bool first_term = true;
  // print nothing if coeff == 0
  if (!almost_zero(m_origin(i))) {
    ss << m_origin(i);
    first_term = false;
  }
 
  for (int j = 0; j < m_to_n.cols(); ++j) {
    double coeff = m_to_n(i, j);
 
    // print nothing if coeff == 0
    if (almost_zero(coeff)) {
      continue;
    }
 
    // if coeff < 0
    if (coeff < 0) {
      if (!first_term) {
        ss << " - " << -coeff << "*" << comp_x_str(j);
      } else {
        ss << coeff << "*" << comp_x_str(j);
      }
    }
 
    // if coeff > 0
    else {
      if (!first_term) {
        ss << " + " << coeff << "*" << comp_x_str(j);
      } else {
        ss << coeff << "*" << comp_x_str(j);
      }
    }
    ss << " ";
 
    first_term = false;
  }
 
  return ss.str();
}
 
std::string CompositionConverter::param_chem_pot_formula(size_type i) const {
  // param_chem_pot = m_to_n.transpose() * chem_pot;
 
  std::stringstream ss;
 
  auto print_chem_pot = [&](int j) {
    return "chem_pot(" + m_components[j] + ") ";
  };
 
  ss << "param_chem_pot(" << comp_var(i) << ") = ";
  Eigen::MatrixXd Mt = m_to_n.transpose();
  bool first_term = true;
  for (int j = 0; j < Mt.cols(); ++j) {
    double coeff = Mt(i, j);
 
    // print nothing if n == 0
 
    if (almost_zero(coeff) || Local::_is_vacancy(m_components[j])) {
      continue;
    }
 
    // print 'A' or '+A' if n == 1
    if (almost_zero(coeff - 1)) {
      if (!first_term) {
        ss << "+ ";
      }
      ss << print_chem_pot(j);
    }
 
    // print '-A' if n == -1
    else if (almost_zero(coeff + 1)) {
      ss << "- " << print_chem_pot(j);
    }
 
    // print 'nA' or '+nA' if n > 0
    else if (coeff > 0) {
      if (!first_term) {
        ss << "+ ";
      }
      ss << coeff << "*" << print_chem_pot(j);
    }
 
    // print '-nA' if n < 0
    else {
      ss << coeff << "*" << print_chem_pot(j);
    }
 
    first_term = false;
  }
 
  return ss.str();
}
 
std::string CompositionConverter::origin_formula() const {
  return _n_formula(m_origin);
}
 
std::string CompositionConverter::end_member_formula(size_type i) const {
  return _n_formula(m_end_members.col(i));
}
 
void CompositionConverter::_add_end_member(Eigen::VectorXd _end_member) {
  _check_size(_end_member);
 
  int r = m_components.size();
  int c = m_end_members.cols();
 
  Eigen::MatrixXd tmp(r, c + 1);
  tmp.leftCols(c) = m_end_members;
  tmp.col(c) = _end_member;
  m_end_members = tmp;
}
 
void CompositionConverter::_check_size(const Eigen::MatrixXd &mat) const {
  if (m_components.size() != mat.rows()) {
    throw std::runtime_error(
        std::string("Error in CompositionConverter: origin or end member "
                    "vector size does not match components size."));
  }
}
 
void CompositionConverter::_calc_conversion_matrices() {
  // calculate m_to_n and m_to_x:
  //
  // n = origin + m_to_n*x
  // x = m_to_x*(n - origin)
 
  // end_members.col(i) corresponds to x such that x[i] = 1, x[j!=i] = 0,
  //  -> end_members.col(i) = origin + m_to_n.col(i)
 
  int c = m_end_members.cols();
 
  m_to_n = m_end_members.leftCols(c).colwise() - m_origin;
 
  // x = m_to_x*(n - origin)
  //   -> x = m_to_x*(origin + m_to_n*x - origin)
  //   -> x = m_to_x*m_to_n*x
 
  // m_to_x is left pseudoinverse of m_to_n, which must have full column rank,
  // because it describes the composition space:
  //
  // I = A+ * A, (A+ is the left pseudoinverse)
  // if A has full column rank, (A.t * A) is invertible, so
  //   A+ = (A.t * A).inv * A.t
  m_to_x = (m_to_n.transpose() * m_to_n).inverse() * m_to_n.transpose();
}
 
std::string CompositionConverter::_n_formula(const Eigen::VectorXd &vec) const {
  std::stringstream tstr;
 
  // for each molecule:
  for (int i = 0; i < vec.size(); i++) {
    // print 'A' if x == 1
    if (almost_zero(vec(i) - 1)) {
      tstr << m_components[i];
    }
    // print 'A(x)' if x != 0
    else if (!almost_zero(vec(i))) {
      tstr << m_components[i] << "(" << vec(i) << ")";
    }
  }
 
  return tstr.str();
}
 
void display_composition_axes(
    std::ostream &stream,
    const std::map<std::string, CompositionConverter> &map) {
  if (map.size() == 0) {
    return;
  }
 
  auto comp_var = CompositionConverter::comp_var;
 
  stream << std::setw(10) << "KEY"
         << " ";
  stream << std::setw(10) << "ORIGIN"
         << " ";
  for (int i = 0; i < map.begin()->second.independent_compositions(); i++) {
    stream << std::setw(10) << comp_var(i) << " ";
  }
  stream << "    ";
  stream << "GENERAL FORMULA";
  stream << std::endl;
 
  stream << std::setw(10) << "  ---"
         << " ";
  stream << std::setw(10) << "  ---"
         << " ";
  for (int i = 0; i < map.begin()->second.independent_compositions(); i++) {
    stream << std::setw(10) << "  ---"
           << " ";
  }
  stream << "    ";
  stream << "---" << std::endl;
 
  for (auto it = map.cbegin(); it != map.cend(); ++it) {
    stream << std::setw(10) << it->first << " ";
    stream << std::setw(10) << it->second.origin_formula() << " ";
    for (int i = 0; i < it->second.independent_compositions(); ++i) {
      stream << std::setw(10) << it->second.end_member_formula(i) << " ";
    }
    stream << "    ";
    stream << std::setw(10) << it->second.mol_formula() << "\n";
  }
}
 
void display_comp(std::ostream &stream, const CompositionConverter &f,
                  int indent) {
  for (int i = 0; i < f.independent_compositions(); ++i) {
    stream << std::string(indent, ' ') << f.comp_formula(i) << "\n";
  }
}
 
void display_comp_n(std::ostream &stream, const CompositionConverter &f,
                    int indent) {
  for (int i = 0; i < f.components().size(); ++i) {
    stream << std::string(indent, ' ') << f.comp_n_formula(i) << "\n";
  }
}
 
void display_param_chem_pot(std::ostream &stream, const CompositionConverter &f,
                            int indent) {
  for (int i = 0; i < f.independent_compositions(); ++i) {
    stream << std::string(indent, ' ') << f.param_chem_pot_formula(i) << "\n";
  }
}
 
Eigen::MatrixXd end_members(
    const ParamComposition::AllowedOccupants &_allowed_occs) {
  ParamComposition param_comp(_allowed_occs);
  param_comp.generate_prim_end_members();
  return param_comp.prim_end_members().transpose();
}
 
Eigen::MatrixXd _composition_space(
    const ParamComposition::AllowedOccupants &_allowed_occs, double tol) {
  // Get Va index if it exists, and store 0 or 1 in N_Va
  std::vector<std::string> struc_mol_name =
      ParamComposition(_allowed_occs).components();
  Index Va_index = find_index_if(struc_mol_name, [=](const std::string &str) {
    return Local::_is_vacancy(str);
  });
  bool has_Va = (Va_index != struc_mol_name.size());
 
  // convert to atom frac
  Eigen::MatrixXd E = end_members(_allowed_occs);
  if (has_Va) {
    E.row(Va_index) = Eigen::VectorXd::Zero(E.cols());
  }
  for (int i = 0; i < E.cols(); i++) {
    E.col(i) /= E.col(i).sum();
  }
 
  // convert to atom frac space
  Eigen::MatrixXd M(E.rows(), E.cols() - 1);
  for (int i = 0; i < M.cols(); ++i) {
    M.col(i) = E.col(i + 1) - E.col(0);
  }
  return M;
}
 
Eigen::MatrixXd composition_space(
    const ParamComposition::AllowedOccupants &_allowed_occs, double tol) {
  auto Qr = _composition_space(_allowed_occs, tol).fullPivHouseholderQr();
  Qr.setThreshold(tol);
  auto Q = Qr.matrixQ();
  return Q.leftCols(Qr.rank());
}
 
Eigen::MatrixXd null_composition_space(
    const ParamComposition::AllowedOccupants &_allowed_occs, double tol) {
  auto Qr = _composition_space(_allowed_occs, tol).fullPivHouseholderQr();
  Qr.setThreshold(tol);
  auto Q = Qr.matrixQ();
  return Q.rightCols(Q.cols() - Qr.rank());
}
 
}  // namespace CASM