OccLocation.cc

Go to the documentation of this file.

#include "casm/monte_carlo/OccLocation.hh"
 
#include "casm/clex/Configuration.hh"
#include "casm/external/MersenneTwister/MersenneTwister.h"
#include "casm/monte_carlo/Conversions.hh"
#include "casm/monte_carlo/OccCandidate.hh"
 
namespace CASM {
namespace Monte {
 
OccLocation::OccLocation(const Conversions &_convert,
                         const OccCandidateList &_cand)
    : m_convert(_convert), m_cand(_cand), m_loc(_cand.size()), m_kmc(false) {}
 
void OccLocation::initialize(const Configuration &config) {
  m_mol.clear();
  m_species.clear();
  m_l_to_mol.clear();
  for (auto &vec : m_loc) {
    vec.clear();
  }
 
  Index Nmut = 0;
  for (Index l = 0; l < config.size(); ++l) {
    Index asym = m_convert.l_to_asym(l);
    if (m_convert.occ_size(asym) > 1) {
      Nmut++;
    }
  }
 
  m_mol.resize(Nmut);
  m_l_to_mol.reserve(config.size());
  Index mol_id = 0;
  for (Index l = 0; l < config.size(); ++l) {
    Index asym = m_convert.l_to_asym(l);
    if (m_convert.occ_size(asym) > 1) {
      Index species_index = m_convert.species_index(asym, config.occ(l));
      Index cand_index = m_cand.index(asym, species_index);
 
      Mol &mol = m_mol[mol_id];
      mol.id = mol_id;
      mol.l = l;
      mol.asym = asym;
      mol.species_index = species_index;
      mol.loc = m_loc[cand_index].size();
 
      if (m_kmc) {
        // only atoms now
        for (Index j = 0; j < 1; j++) {
          Species spec(m_convert.l_to_bijk(l));
          spec.species_index = species_index;
          spec.id = m_species.size();
          mol.component.push_back(spec.id);
 
          m_species.push_back(spec);
        }
      }
 
      m_loc[cand_index].push_back(mol_id);
      m_l_to_mol.push_back(mol_id);
      mol_id++;
    } else {
      m_l_to_mol.push_back(Nmut);
    }
  }
  if (m_kmc) {
    m_tmol = m_mol;
  }
}
 
OccLocation::size_type OccLocation::size() const { return m_mol.size(); }
 
Mol &OccLocation::mol(Index mol_id) { return m_mol[mol_id]; }
 
const Mol &OccLocation::mol(Index mol_id) const { return m_mol[mol_id]; }
 
OccLocation::size_type OccLocation::cand_size(Index cand_index) const {
  return m_loc[cand_index].size();
}
 
OccLocation::size_type OccLocation::cand_size(const OccCandidate &cand) const {
  return cand_size(m_cand.index(cand));
}
 
Index OccLocation::mol_id(Index cand_index, Index loc) const {
  return m_loc[cand_index][loc];
}
 
Index OccLocation::mol_id(const OccCandidate &cand, Index loc) const {
  return mol_id(m_cand.index(cand), loc);
}
 
Index OccLocation::l_to_mol_id(Index l) const { return m_l_to_mol[l]; }
 
OccEvent &OccLocation::propose_canonical(
    OccEvent &e, const std::vector<OccSwap> &canonical_swap,
    MTRand &mtrand) const {
  Index tsize = canonical_swap.size();
  m_tsum.resize(tsize + 1);
 
  m_tsum[0] = 0.;
  for (Index i = 0; i < tsize; ++i) {
    m_tsum[i + 1] =
        m_tsum[i] + ((double)cand_size(canonical_swap[i].cand_a)) *
                        ((double)cand_size(canonical_swap[i].cand_b));
  }
 
  double rand = mtrand.randExc(m_tsum.back());
 
  for (Index i = 0; i < tsize; ++i) {
    if (rand < m_tsum[i + 1]) {
      return _propose(e, canonical_swap[i], mtrand);
    }
  }
 
  throw std::runtime_error("OccLocation::propose_canonical error");
}
 
OccEvent &OccLocation::propose_grand_canonical(OccEvent &e, const OccSwap &swap,
                                               MTRand &mtrand) const {
  e.occ_transform.resize(1);
  e.species_traj.resize(0);
 
  Index index_cand_a = m_cand.index(swap.cand_a);
  Index index_cand_b = m_cand.index(swap.cand_b);
 
  OccTransform &f_a = e.occ_transform[0];
  f_a.mol_id = m_loc[index_cand_a][mtrand.randInt(cand_size(swap.cand_a) - 1)];
  f_a.l = m_mol[f_a.mol_id].l;
  f_a.asym = m_cand[index_cand_a].asym;
  f_a.from_species = m_cand[index_cand_a].species_index;
  f_a.to_species = m_cand[index_cand_b].species_index;
 
  return e;
}
 
void OccLocation::apply(const OccEvent &e, ConfigDoF &configdof) {
  // copy original Mol.component
  if (m_kmc) {
    for (const auto &occ : e.occ_transform) {
      m_tmol[occ.mol_id].component = m_mol[occ.mol_id].component;
    }
  }
 
  // update Mol and config occupation
  for (const auto &occ : e.occ_transform) {
    auto &mol = m_mol[occ.mol_id];
 
    // set config occupation
    configdof.occ(mol.l) = m_convert.occ_index(mol.asym, occ.to_species);
 
    // remove from m_loc
    Index cand_index = m_cand.index(mol.asym, mol.species_index);
    Index back = m_loc[cand_index].back();
    m_loc[cand_index][mol.loc] = back;
    m_mol[back].loc = mol.loc;
    m_loc[cand_index].pop_back();
 
    // set Mol.species index
    mol.species_index = occ.to_species;
 
    if (m_kmc) {
      mol.component.resize(m_convert.components_size(mol.species_index));
    }
 
    // add to m_loc
    cand_index = m_cand.index(mol.asym, mol.species_index);
    mol.loc = m_loc[cand_index].size();
    m_loc[cand_index].push_back(mol.id);
  }
 
  if (m_kmc) {
    // update Mol.component
    for (const auto &traj : e.species_traj) {
      m_mol[traj.to.mol_id].component[traj.to.mol_comp] =
          m_tmol[traj.from.mol_id].component[traj.from.mol_comp];
    }
  }
}
 
OccEvent &OccLocation::_propose(OccEvent &e, const OccSwap &swap,
                                MTRand &mtrand, Index cand_a, Index cand_b,
                                Index size_a, Index size_b) const {
  e.occ_transform.resize(2);
  e.species_traj.resize(0);
 
  OccTransform &f_a = e.occ_transform[0];
  f_a.mol_id = m_loc[cand_a][mtrand.randInt(size_a - 1)];
  f_a.l = m_mol[f_a.mol_id].l;
  f_a.asym = m_cand[cand_a].asym;
  f_a.from_species = m_cand[cand_a].species_index;
  f_a.to_species = m_cand[cand_b].species_index;
  // std::cout << "size_a: " << size_a << "  loc: " << m_mol[f_a.mol_id].loc <<
  // std::endl;
 
  OccTransform &f_b = e.occ_transform[1];
  f_b.mol_id = m_loc[cand_b][mtrand.randInt(size_b - 1)];
  f_b.l = m_mol[f_b.mol_id].l;
  f_b.asym = m_cand[cand_b].asym;
  f_b.from_species = m_cand[cand_b].species_index;
  f_b.to_species = m_cand[cand_a].species_index;
  // std::cout << "size_b: " << size_b << "  loc: " << m_mol[f_b.mol_id].loc <<
  // std::endl;
 
  return e;
}
 
OccEvent &OccLocation::_propose(OccEvent &e, const OccSwap &swap,
                                MTRand &mtrand) const {
  Index cand_a = m_cand.index(swap.cand_a);
  Index cand_b = m_cand.index(swap.cand_b);
  Index size_a = m_loc[cand_a].size();
  Index size_b = m_loc[cand_b].size();
  return _propose(e, swap, mtrand, cand_a, cand_b, size_a, size_b);
}
 
}  // namespace Monte
}  // namespace CASM