SiteExchanger.cc

Go to the documentation of this file.

#include "casm/monte_carlo/SiteExchanger.hh"
 
#include "casm/casm_io/container/stream_io.hh"
#include "casm/clex/CompositionConverter.hh"
#include "casm/clex/PrimClex.hh"
#include "casm/clex/Supercell.hh"
#include "casm/crystallography/Structure.hh"
 
namespace CASM {
namespace Monte {
 
// ---- SiteExchanger Definitions ---------------------------------
 
SiteExchanger::SiteExchanger(const Supercell &scel) {
  // The occupants in a Configuration are ordered in blocks of basis sites.
  int scel_volume = scel.volume();
 
  // int scel_basis = scel->basis_size();
  int prim_basis = scel.prim().basis().size();
  std::vector<std::string> allowed_components =
      scel.primclex().composition_axes().components();
 
  // Count over sites in prim basis.
  for (Index prim_basis_site = 0; prim_basis_site < prim_basis;
       prim_basis_site++) {
    // If the site we're working allows multiple occupants, we're interested in
    // filling up values for it.
    const auto &allowed =
        scel.prim().basis()[prim_basis_site].allowed_occupants();
    std::vector<std::string> site_allowed_occ(allowed.cbegin(), allowed.cend());
 
    if (site_allowed_occ.size() > 1) {
      // This is the center array of possible_swap and changes with each
      // prim_basis_site with more than one allowed occupant
      std::vector<std::vector<int> > single_possible_swap;
 
      // This is a site_to_mol inner array and changes with each prim_basis_site
      // with more than one allowed occupant
      std::vector<int> single_site_to_mol;
 
      // For every current occupant we could have, we need an std::vector<int>
      // for single_possible_swap and an int for single_site_to_mol
      for (int possible_curr_occ = 0;
           possible_curr_occ < site_allowed_occ.size(); possible_curr_occ++) {
        // This is the inner array of possible_swap and tells you all the OTHER
        // values the current occupant could have been for the site
        std::vector<int> swaps;
 
        // Loop over the same thing, and only keep different values of the
        // current occupant as swaps that could occur
        for (int other_curr_occ = 0; other_curr_occ < site_allowed_occ.size();
             other_curr_occ++) {
          if (possible_curr_occ != other_curr_occ) {
            swaps.push_back(other_curr_occ);
          }
        }
 
        // At this point we have all the available swaps for one site of a given
        // current occupant, so we save them
        single_possible_swap.push_back(swaps);
 
        // For site_to_mol we want a way to translate what possible_curr_occ
        // means in terms of allowed components at the PrimClex level
        int mol_ind =
            std::find(allowed_components.cbegin(), allowed_components.cend(),
                      site_allowed_occ[possible_curr_occ]) -
            allowed_components.cbegin();
 
        // This should never ever happen
        if (mol_ind == allowed_components.size()) {
          std::cerr << "ERROR in Monte::populate_occ_exchange_tables"
                    << std::endl;
          std::cerr << "The possible components are " << allowed_components
                    << std::endl;
          std::cerr << "Could not find "
                    << allowed_components[possible_curr_occ]
                    << " in the allowed components of the PrimClex prim."
                    << std::endl;
          exit(9000);
        }
 
        // Save the index into allowed_components
        single_site_to_mol.push_back(mol_ind);
      }
 
      // possible_swap has blocks of identical double arrays
      // single_possible_swap. They get pushed back the appropriate amount of
      // times here
      m_possible_swap.push_back(single_possible_swap);
 
      // site_to_mol has blocks of identical single arrays single_site_to_mol.
      // They get pushed back the appropriate amount of times here
      m_sublat_to_mol.push_back(single_site_to_mol);
 
      // This loop should happen scel_volume times and takes care of repeating
      // values in m_variable_sites and m_sublat the appropriate amount of
      // times. Both m_variable_sites and m_sublat have the same length on the
      // outside: one slot for each site in the Configuration that can hold more
      // than one occupant.
      for (int variable_site = prim_basis_site * scel_volume;
           variable_site < prim_basis_site * scel_volume + scel_volume;
           variable_site++) {
        // variable_sites determined by the counter. Contains indexes of sites
        // in the Configuration that allow more than one occupant
        m_variable_sites.push_back(variable_site);
 
        // store the sublat index
        m_sublat.push_back(prim_basis_site);
      }
    }
 
    // Here we populate m_possible_swap and m_sublat_to_mol for sublattices that
    // only have a single allowed occupant.
    else {
      // In this else block, single_site_to_mol can only have a single value in
      // the array
      std::vector<int> single_site_to_mol;
 
      // The one and only allowed occupant at prim_basis_site must be
      std::string only_site_occ =
          scel.prim().basis()[prim_basis_site].allowed_occupants()[0];
 
      // And the corresponding index for that occupant in terms of
      // allowed_components is
      int mol_ind = std::find(allowed_components.cbegin(),
                              allowed_components.cend(), only_site_occ) -
                    allowed_components.cbegin();
      if (mol_ind == allowed_components.size()) {
        std::cerr << "ERROR in Monte::populate_occ_exchange_tables"
                  << std::endl;
        std::cerr << "The possible components are " << allowed_components
                  << std::endl;
        std::cerr << "Could not find " << only_site_occ
                  << " in the allowed components of the PrimClex prim."
                  << std::endl;
        exit(9000);
      }
 
      // No swaps are possible on this sublattice, so we push back empty
      // std::vector<std::vector<int> >
      std::vector<std::vector<int> > single_possible_swap;
      m_possible_swap.push_back(single_possible_swap);
 
      // Put that single value we found into the array
      single_site_to_mol.push_back(mol_ind);
 
      m_sublat_to_mol.push_back(single_site_to_mol);
    }
  }
 
  return;
}
 
}  // namespace Monte
}  // namespace CASM