GrandCanonicalIO.cc

Go to the documentation of this file.

#include "casm/monte_carlo/grand_canonical/GrandCanonicalIO.hh"
 
#include "casm/casm_io/dataformatter/DataFormatter.hh"
#include "casm/crystallography/BasicStructure.hh"
#include "casm/global/definitions.hh"
#include "casm/monte_carlo/MonteIO_impl.hh"
#include "casm/monte_carlo/grand_canonical/GrandCanonical.hh"
 
namespace CASM {
namespace Monte {
 
DataFormatter<ConstMonteCarloPtr> make_results_formatter(
    const GrandCanonical &mc) {
  DataFormatter<ConstMonteCarloPtr> formatter;
 
  formatter.push_back(MonteCarloIsEquilibratedFormatter());
  formatter.push_back(MonteCarloNEquilSamplesFormatter());
  formatter.push_back(MonteCarloIsConvergedFormatter());
  formatter.push_back(MonteCarloNAvgSamplesFormatter());
 
  formatter.push_back(MonteCarloTFormatter<GrandCanonical>());
 
  std::set<std::string> exclude;
  std::string name;
 
  // always sample Beta, potential_energy, and formation_energy
  {
    formatter.push_back(MonteCarloBetaFormatter<GrandCanonical>());
    name = "potential_energy";
    formatter.push_back(MonteCarloMeanFormatter(name));
    formatter.push_back(MonteCarloPrecFormatter(name));
    exclude.insert(name);
 
    name = "formation_energy";
    formatter.push_back(MonteCarloMeanFormatter(name));
    formatter.push_back(MonteCarloPrecFormatter(name));
    exclude.insert(name);
  }
 
  // always sample param_chem_pot, comp
  for (int i = 0;
       i < mc.primclex().composition_axes().independent_compositions(); ++i) {
    formatter.push_back(MonteCarloParamChemPotFormatter<GrandCanonical>(mc, i));
  }
 
  for (int i = 0;
       i < mc.primclex().composition_axes().independent_compositions(); i++) {
    name = std::string("comp(") + mc.primclex().composition_axes().comp_var(i) +
           ")";
    formatter.push_back(MonteCarloMeanFormatter(name));
    formatter.push_back(MonteCarloPrecFormatter(name));
    exclude.insert(name);
  }
 
  // always sample comp_n
  auto struc_mol_name = xtal::struc_molecule_name(mc.primclex().prim());
  for (int i = 0; i < struc_mol_name.size(); ++i) {
    name = std::string("comp_n(") + struc_mol_name[i] + ")";
    formatter.push_back(MonteCarloMeanFormatter(name));
    formatter.push_back(MonteCarloPrecFormatter(name));
    exclude.insert(name);
  }
 
  // include mean/prec of other properties
  for (auto it = mc.samplers().cbegin(); it != mc.samplers().cend(); ++it) {
    if (exclude.find(it->first) == exclude.end()) {
      formatter.push_back(MonteCarloMeanFormatter(it->first));
      formatter.push_back(MonteCarloPrecFormatter(it->first));
    }
  }
 
  // include heat_capacity
  formatter.push_back(MonteCarloHeatCapacityFormatter<GrandCanonical>());
 
  // include susc_x
  for (int i = 0;
       i < mc.primclex().composition_axes().independent_compositions(); i++) {
    for (int j = i;
         j < mc.primclex().composition_axes().independent_compositions(); j++) {
      auto comp_var_i = mc.primclex().composition_axes().comp_var(i);
      auto comp_var_j = mc.primclex().composition_axes().comp_var(j);
      formatter.push_back(
          MonteCarloSuscXFormatter<GrandCanonical>(comp_var_i, comp_var_j));
    }
  }
 
  // include thermo-chem susc
  for (int i = 0;
       i < mc.primclex().composition_axes().independent_compositions(); i++) {
    auto comp_var_i = mc.primclex().composition_axes().comp_var(i);
    formatter.push_back(
        MonteCarloThermoChemSuscXFormatter<GrandCanonical>(comp_var_i));
  }
 
  // include susc_n
  for (int i = 0; i < struc_mol_name.size(); ++i) {
    for (int j = i; j < struc_mol_name.size(); ++j) {
      auto species_i = struc_mol_name[i];
      auto species_j = struc_mol_name[j];
      formatter.push_back(
          MonteCarloSuscNFormatter<GrandCanonical>(species_i, species_j));
    }
  }
 
  // include thermo-chem susc
  for (int i = 0; i < struc_mol_name.size(); ++i) {
    auto species_i = struc_mol_name[i];
    formatter.push_back(
        MonteCarloThermoChemSuscNFormatter<GrandCanonical>(species_i));
  }
 
  return formatter;
}
 
DataFormatter<ConstMonteCarloPtr> make_lte_results_formatter(
    const GrandCanonical &mc, const double &phi_LTE1,
    const std::string &configname) {
  DataFormatter<ConstMonteCarloPtr> formatter;
 
  bool print_json = true;
  formatter.push_back(ConstantValueFormatter<std::string, ConstMonteCarloPtr>(
      "configname", configname, print_json));
  formatter.push_back(MonteCarloTFormatter<GrandCanonical>());
  formatter.push_back(GrandCanonicalLTEFormatter(phi_LTE1));
  std::set<std::string> exclude;
  std::string name;
 
  // always sample Beta, potential_energy, and formation_energy
  {
    formatter.push_back(MonteCarloBetaFormatter<GrandCanonical>());
    name = "gs_potential_energy";
    auto evaluator = [=](const ConstMonteCarloPtr &ptr) {
      return static_cast<const GrandCanonical *>(ptr)->potential_energy();
    };
    formatter.push_back(GenericDatumFormatter<double, ConstMonteCarloPtr>(
        name, name, evaluator));
    exclude.insert(name);
  }
 
  {
    name = "gs_formation_energy";
    auto evaluator = [=](const ConstMonteCarloPtr &ptr) {
      return static_cast<const GrandCanonical *>(ptr)->formation_energy();
    };
    formatter.push_back(GenericDatumFormatter<double, ConstMonteCarloPtr>(
        name, name, evaluator));
    exclude.insert(name);
  }
 
  // always sample param_chem_pot, comp
  for (int i = 0;
       i < mc.primclex().composition_axes().independent_compositions(); ++i) {
    formatter.push_back(MonteCarloParamChemPotFormatter<GrandCanonical>(mc, i));
  }
 
  for (int i = 0;
       i < mc.primclex().composition_axes().independent_compositions(); i++) {
    name = std::string("gs_comp(") +
           mc.primclex().composition_axes().comp_var(i) + ")";
    auto evaluator = [=](const ConstMonteCarloPtr &ptr) {
      const GrandCanonical *_ptr = static_cast<const GrandCanonical *>(ptr);
      return _ptr->primclex().composition_axes().param_composition(
          _ptr->comp_n())[i];
    };
    formatter.push_back(GenericDatumFormatter<double, ConstMonteCarloPtr>(
        name, name, evaluator));
    exclude.insert(name);
  }
 
  // always sample comp_n
  auto struc_mol_name = xtal::struc_molecule_name(mc.primclex().prim());
  for (int i = 0; i < struc_mol_name.size(); ++i) {
    name = std::string("gs_comp_n(") + struc_mol_name[i] + ")";
    auto evaluator = [=](const ConstMonteCarloPtr &ptr) {
      return static_cast<const GrandCanonical *>(ptr)->comp_n()[i];
    };
    formatter.push_back(GenericDatumFormatter<double, ConstMonteCarloPtr>(
        name, name, evaluator));
    exclude.insert(name);
  }
 
  return formatter;
}
 
jsonParser &to_json(const GrandCanonicalConditions &conditions,
                    jsonParser &json) {
  json.put_obj();
  json["temperature"] = conditions.temperature();
  json["param_chem_pot"] = jsonParser::object();
  auto param_chem_pot = conditions.param_chem_pot();
  for (int i = 0; i < param_chem_pot.size(); i++) {
    json["param_chem_pot"][CompositionConverter::comp_var(i)] =
        param_chem_pot(i);
  }
  json["tolerance"] = conditions.tolerance();
 
  return json;
}
 
void from_json(GrandCanonicalConditions &conditions, const PrimClex &primclex,
               const jsonParser &json) {
  double temp = json["temperature"].get<double>();
  double tol = json["tolerance"].get<double>();
 
  int Nparam = primclex.composition_axes().independent_compositions();
  Eigen::VectorXd param_chem_pot(Nparam);
 
  for (int i = 0; i < Nparam; i++) {
    param_chem_pot(i) =
        json["param_chem_pot"][CompositionConverter::comp_var(i)].get<double>();
  }
 
  conditions = GrandCanonicalConditions(primclex, temp, param_chem_pot, tol);
}
 
GenericDatumFormatter<double, ConstMonteCarloPtr> GrandCanonicalLTEFormatter(
    const double &phi_LTE1) {
  auto evaluator = [=](const ConstMonteCarloPtr &mc) { return phi_LTE1; };
  return GenericDatumFormatter<double, ConstMonteCarloPtr>("phi_LTE", "phi_LTE",
                                                           evaluator);
}
 
void write_lte_results(const MonteSettings &settings, const GrandCanonical &mc,
                       const double &phi_LTE1, const std::string &configname,
                       Log &_log) {
  try {
    fs::create_directories(settings.output_directory());
    MonteCarloDirectoryStructure dir(settings.output_directory());
    auto formatter = make_lte_results_formatter(mc, phi_LTE1, configname);
 
    // write csv path results
    if (settings.write_csv()) {
      fs::path file = dir.results_csv();
      _log << "write: " << dir.results_csv() << "\n";
      fs::ofstream sout;
 
      if (!fs::exists(file)) {
        sout.open(file);
        formatter.print_header(&mc, sout);
      } else {
        sout.open(file, std::ios::app);
      }
 
      formatter.print(&mc, sout);
 
      sout.close();
    }
 
    // write json path results
    if (settings.write_json()) {
      fs::path file = dir.results_json();
      _log << "write: " << dir.results_json() << "\n";
 
      jsonParser results;
      if (fs::exists(file)) {
        results.read(file);
      } else {
        results = jsonParser::object();
      }
 
      formatter.to_json_arrays(&mc, results);
      results.write(file);
    }
  } catch (...) {
    std::cerr << "ERROR writing results" << std::endl;
    throw;
  }
}
 
}  // namespace Monte
}  // namespace CASM