GrandCanonicalConditions.cc

Go to the documentation of this file.

#include "casm/monte_carlo/grand_canonical/GrandCanonicalConditions.hh"
 
#include "casm/casm_io/container/json_io.hh"
#include "casm/clex/PrimClex.hh"
#include "casm/monte_carlo/MonteCarlo.hh"
#include "casm/monte_carlo/MonteSettings.hh"
 
namespace CASM {
namespace Monte {
 
GrandCanonicalConditions::GrandCanonicalConditions(
    const PrimClex &_primclex, double _temperature,
    const Eigen::VectorXd &_param_chem_pot, double _tol)
    :
 
      m_primclex(&_primclex),
      m_tolerance(_tol) {
  // -- set T ----
  set_temperature(_temperature);
 
  // -- set param_chem_pot ----
  set_param_chem_pot(_param_chem_pot);
}
 
// ***************************************ACCESSORS**********************************************
// //
 
const PrimClex &GrandCanonicalConditions::primclex() const {
  return *m_primclex;
}
 
double GrandCanonicalConditions::temperature() const { return m_temperature; }
 
double GrandCanonicalConditions::beta() const { return m_beta; }
 
Eigen::MatrixXd GrandCanonicalConditions::exchange_chem_pot() const {
  return m_exchange_chem_pot;
}
 
double GrandCanonicalConditions::exchange_chem_pot(Index index_new,
                                                   Index index_curr) const {
  return m_exchange_chem_pot(index_new, index_curr);
}
 
Eigen::VectorXd GrandCanonicalConditions::param_chem_pot() const {
  return m_param_chem_pot;
}
 
double GrandCanonicalConditions::param_chem_pot(Index index) const {
  return m_param_chem_pot(index);
}
 
double GrandCanonicalConditions::tolerance() const { return m_tolerance; }
 
// ***************************************MUTATORS***********************************************
// //
 
void GrandCanonicalConditions::set_temperature(double in_temp) {
  m_temperature = in_temp;
  m_beta = 1.0 / (KB * m_temperature);
  return;
}
 
void GrandCanonicalConditions::set_param_chem_pot(
    const Eigen::VectorXd &in_param_chem_pot) {
  m_param_chem_pot = in_param_chem_pot;
  m_chem_pot = primclex().composition_axes().dparam_dmol().transpose() *
               m_param_chem_pot;
 
  int Ncomp = primclex().composition_axes().components().size();
  m_exchange_chem_pot = Eigen::MatrixXd(Ncomp, Ncomp);
  for (int index_new = 0; index_new < Ncomp; ++index_new) {
    for (int index_curr = 0; index_curr < Ncomp; ++index_curr) {
      Eigen::VectorXl dn = Eigen::VectorXl::Zero(Ncomp);
      dn(index_new) += 1;
      dn(index_curr) -= 1;
      m_exchange_chem_pot(index_new, index_curr) =
          m_param_chem_pot.transpose() *
          primclex().composition_axes().dparam_dmol() * dn.cast<double>();
    }
  }
 
  return;
}
 
void GrandCanonicalConditions::set_param_chem_pot(Index ind,
                                                  double in_param_chem_pot) {
  m_param_chem_pot(ind) = in_param_chem_pot;
  set_param_chem_pot(m_param_chem_pot);
  return;
}
 
// ***************************************OPERATORS**********************************************
// //
 
GrandCanonicalConditions &GrandCanonicalConditions::operator+=(
    const GrandCanonicalConditions &RHS) {
  m_temperature += RHS.m_temperature;
  for (int i = 0; i < m_param_chem_pot.size(); i++) {
    m_param_chem_pot(i) += RHS.m_param_chem_pot(i);
  }
  set_param_chem_pot(m_param_chem_pot);
  m_beta = 1.0 / (CASM::KB * m_temperature);
  return *this;
}
 
GrandCanonicalConditions GrandCanonicalConditions::operator+(
    const GrandCanonicalConditions &RHS) const {
  return GrandCanonicalConditions(*this) += RHS;
}
 
GrandCanonicalConditions &GrandCanonicalConditions::operator-=(
    const GrandCanonicalConditions &RHS) {
  m_temperature -= RHS.m_temperature;
  for (int i = 0; i < m_param_chem_pot.size(); i++) {
    m_param_chem_pot(i) -= RHS.m_param_chem_pot(i);
  }
  set_param_chem_pot(m_param_chem_pot);
  m_beta = 1.0 / (CASM::KB * m_temperature);
  return *this;
}
 
GrandCanonicalConditions GrandCanonicalConditions::operator-(
    const GrandCanonicalConditions &RHS) const {
  return GrandCanonicalConditions(*this) -= RHS;
}
 
bool GrandCanonicalConditions::operator==(
    const GrandCanonicalConditions &RHS) const {
  if (!almost_zero(m_temperature - RHS.m_temperature, m_tolerance)) {
    return false;
  }
 
  for (int i = 0; i < m_param_chem_pot.size(); i++) {
    if (!almost_zero(m_param_chem_pot(i) - RHS.m_param_chem_pot(i),
                     m_tolerance)) {
      return false;
    }
  }
 
  return true;
}
 
bool GrandCanonicalConditions::operator!=(
    const GrandCanonicalConditions &RHS) const {
  return !(*this == RHS);
}
 
int GrandCanonicalConditions::operator/(
    const GrandCanonicalConditions &RHS_inc) const {
  int max_division = 0;
 
  if (!almost_zero(RHS_inc.temperature())) {
    max_division = round(temperature() / RHS_inc.temperature());
  }
 
  for (Index i = 0; i < param_chem_pot().size(); i++) {
    int temp_division = round(param_chem_pot(i) / RHS_inc.param_chem_pot(i));
 
    if (temp_division > max_division &&
        !almost_zero(RHS_inc.param_chem_pot(i))) {
      max_division = temp_division;
    }
  }
 
  return max_division;
}
 
std::ostream &operator<<(std::ostream &sout,
                         const GrandCanonicalConditions &cond) {
  sout << "T: " << cond.temperature() << "\n";
  for (int i = 0; i < cond.param_chem_pot().size(); i++) {
    jsonParser json;
    sout << "param_chem_pot: " << to_json_array(cond.param_chem_pot(), json)
         << "\n";
  }
  return sout;
}
 
}  // namespace Monte
}  // namespace CASM