MonteCarlo.hh

Go to the documentation of this file.

#ifndef CASM_MonteCarlo_HH
#define CASM_MonteCarlo_HH

#include <vector>
#include "casm/misc/cloneable_ptr.hh"
#include "casm/external/MersenneTwister/MersenneTwister.h"
#include "casm/casm_io/Log.hh"
#include "casm/clex/PrimClex.hh"
#include "casm/clex/Supercell.hh"
#include "casm/clex/ConfigDoF.hh"
#include "casm/monte_carlo/MonteSettings.hh"
#include "casm/monte_carlo/MonteSampler.hh"
#include "casm/monte_carlo/MonteCounter.hh"

namespace CASM {

  struct SamplerNameCompare {

    SamplerNameCompare() {};

    bool operator()(const std::string &A, const std::string &B) const;

  };


  class MonteCarlo {


  public:

    typedef MonteSampler::size_type size_type;

    typedef std::map<std::string, double> ScalarPropertyMap;

    typedef std::map<std::string, Eigen::VectorXd> VectorPropertyMap;

    typedef std::map<std::string, notstd::cloneable_ptr<MonteSampler>, SamplerNameCompare> SamplerMap;

    typedef std::vector<std::pair<MonteCounter::size_type, MonteCounter::size_type> > SampleTimes;


    // ---- Accessors -----------------------------

    const MonteSettings &settings() const {
      return m_settings;
    }

    const PrimClex &primclex() const {
      return m_primclex;
    }

    const Supercell &supercell() const {
      return m_scel;
    }

    void set_nlist() {
      m_nlist = &(supercell().nlist());
    }

    const SuperNeighborList &nlist() const {
      return *m_nlist;
    }

    const Configuration &config() const {
      return m_config;
    }

    const ConfigDoF &configdof() const {
      return m_configdof;
    }


    // ---- Accessors -----------------------------

    void reset(const ConfigDoF &dof) {
      _configdof() = dof;
      clear_samples();
    }


    // ---- Properties ----------------

    const ScalarPropertyMap &scalar_properties() const {
      return m_scalar_property;
    }

    const double &scalar_property(std::string property_name) const {
      return m_scalar_property.find(property_name)->second;
    }

    const VectorPropertyMap &vector_properties() const {
      return m_vector_property;
    }

    const Eigen::VectorXd &vector_property(std::string property_name) const {
      return m_vector_property.find(property_name)->second;
    }


    // ---- Data sampling -------------

    void sample_data(const MonteCounter &counter);

    void clear_samples();

    bool must_converge() const {
      return m_must_converge;
    }

    std::pair<bool, MonteSampler::size_type> is_equilibrated() const;

    bool check_convergence_time() const;

    bool is_converged() const;

    const SamplerMap &samplers() const {
      return m_sampler;
    }

    const SampleTimes &sample_times() const {
      return m_sample_time;
    }

    const std::vector<ConfigDoF> &trajectory() const {
      return m_trajectory;
    }

    bool debug() const {
      return m_debug;
    }


  protected:

    template<typename MonteTypeSettings>
    MonteCarlo(PrimClex &primclex, const MonteTypeSettings &settings, Log &_log);

    PrimClex &_primclex() const {
      return m_primclex;
    }

    Supercell &_supercell() const {
      return m_scel;
    }

    Configuration &_config() const {
      return m_config;
    }

    ConfigDoF &_configdof() const {
      return m_configdof;
    }

    Log &_log() const {
      return m_log;
    }

    MTRand &_mtrand() {
      return m_twister;
    }

    ScalarPropertyMap &_scalar_properties() {
      return m_scalar_property;
    }

    double &_scalar_property(std::string property_name) {
      return m_scalar_property.find(property_name)->second;
    }

    VectorPropertyMap &_vector_properties() {
      return m_vector_property;
    }

    Eigen::VectorXd &_vector_property(std::string property_name) {
      return m_vector_property.find(property_name)->second;
    }

  private:

    ScalarPropertyMap m_scalar_property;


    VectorPropertyMap m_vector_property;

    const MonteSettings &m_settings;

    PrimClex &m_primclex;

    mutable Supercell m_scel;

    const SuperNeighborList *m_nlist;

    mutable Configuration m_config;

    ConfigDoF &m_configdof;

    MTRand m_twister;


    bool m_write_trajectory = false;

    std::vector<ConfigDoF> m_trajectory;

    bool m_must_converge;

    Log &m_log;

    void _set_check_convergence_time() const;

    SamplerMap m_sampler;

    SampleTimes m_sample_time;


    // Members to remember results of is_equilibrated and is_converged

    mutable std::pair<bool, MonteSampler::size_type> m_is_equil;
    mutable bool m_is_equil_uptodate = false;
    mutable bool m_is_converged;
    mutable bool m_is_converged_uptodate = false;

    // which sample to check convergence after
    mutable MonteSampler::size_type m_next_convergence_check = 100;
    MonteSampler::size_type m_convergence_check_period = 100;

    // in debug mode, allow printing or checking extra things
    bool m_debug;

  };

  template<typename MonteTypeSettings>
  MonteCarlo::MonteCarlo(PrimClex &primclex, const MonteTypeSettings &settings, Log &_log) :
    m_settings(settings),
    m_primclex(primclex),
    m_scel(&primclex, settings.simulation_cell_matrix()),
    m_config(m_scel),
    m_configdof(m_config.configdof()),
    m_write_trajectory(settings.write_trajectory()),
    m_log(_log),
    m_debug(m_settings.debug()) {

    settings.samplers(primclex, std::inserter(m_sampler, m_sampler.begin()));

    m_must_converge = false;
    for(auto it = m_sampler.cbegin(); it != m_sampler.cend(); ++it) {
      if(it->second->must_converge()) {
        m_must_converge = true;
        break;
      }
    }
  }

}
#endif