Developers Guide

Developer's Guide

ここでは、新たに分子力場を追加するための説明をします。

新規な結合長ポテンシャル

結合長にかかるポテンシャルとして、新規な関数を導入したいとしましょう。例えば、結合距離rに対して以下のような関数を使ってみることにします。

V(r) = k * (r - r_0)^4

結合長に関する相互作用はすでに実装されているので（詳細はDesignを参照）、必要なのはこの関数系を表現するPotentialクラスだけです。そしてそのクラスは、この関数のパラメータを格納し、関数の値とその微分の値を計算できる必要があります。

ポテンシャルの実装

クラスを作るに当たって、異なる実数値型のために複数回同じコードを書かなくて済むように、templateパラメータを受け取る準備をしましょう。Mjolnirでは、それらのパラメータはTraitsクラスとして一括で渡されるため、それを一つ受け取って、その中で定義されている型のエイリアスを宣言しておきます。

template<typename traitsT>
class QuarticPotential
{
  public:
    typedef typename traitsT::real_type real_type;
};

エイリアスを宣言する主な理由は、今後メンバ関数やメンバ変数の定義でreal_typeを使う時にいちいちtraitsT::と書きたくないからです。

この関数はパラメータを2つ取ります。kとr_0です。なので、それを持たせましょう。それらを初期化するコンストラクタも作っておきましょう。

template<typename traitsT>
class QuarticPotential
{
  public:
    typedef typename traitsT::real_type real_type;

    QuarticPotential(const real_type k, const real_type r_0)
        : k_(k), r_0_(r_0)
    {}

  private:
    real_type k_, r_0_;
};

では、いよいよ本題に入ります。この関数の値を計算させましょう。あと、微分もできるようにする必要があります。

これらの関数はポテンシャルのパラメータを変えたりしませんから、constをつけておきます。また、例外を投げることもないので、noexceptもつけておきます。

template<typename traitsT>
class QuarticPotential
{
  public:
    typedef typename traitsT::real_type real_type;

    QuarticPotential(const real_type k, const real_type r_0)
        : k_(k), r_0_(r_0)
    {}

    // 変数 r に対して、k * (r-r0)^4 の値を計算する
    real_type potential(const real_type r) const noexcept
    {
        const real_type dr = r - r_0_;
        const real_type dr2 = dr * dr;
        return k_ * dr2 * dr2;
    }

    // 変数 r に対して、ポテンシャルの微分 4k * (r-r0)^3 の値を計算する
    real_type derivative(const real_type r) const noexcept
    {
        const real_type dr = r - r_0_;
        return 4 * k_ * dr * dr * dr;
    }

  private:
    real_type k_, r_0_;
};

残りは、系全体のパラメータ（温度やイオン濃度など）が変化した時のためのupdate関数と、このポテンシャルの名前を取得するためのname関数だけです。今回はこれらの関数は特に何もしないので、そのまま書いてしまいましょう。

template<typename traitsT>
class QuarticPotential
{
  public:
    typedef typename traitsT::real_type real_type;

    QuarticPotential(const real_type k, const real_type r_0)
        : k_(k), r_0_(r_0)
    {}

    // 変数 r に対して、k * (r-r0)^4 の値を計算する
    real_type potential(const real_type r) const noexcept
    {
        const real_type dr = r - r_0_;
        const real_type dr2 = dr * dr;
        return k_ * dr2 * dr2;
    }

    // 変数 r に対して、ポテンシャルの微分 4k * (r-r0)^3 の値を計算する
    real_type derivative(const real_type r) const noexcept
    {
        const real_type dr = r - r_0_;
        return 4 * k_ * dr * dr * dr;
    }

    void update(const system_type&, const real_type) const noexcept {return;}
    const char* name() const noexcept {return "Quartic";}

  private:
    real_type k_, r_0_;
};

このようなクラスの実例には、例えばmjolnir/potential/local/HarmonicPotential.hppを見てください。

ファイル入力

次は、入力ファイルからこのポテンシャルを指定できるようにしましょう。 そのためには、mjolnir/input/read_potential.hppと、mjolnir/input/read_interaction.hppを編集する必要があります。

read_potential.hppでは、それぞれのポテンシャルのパラメータを読むための方法が定義されています。read_interaction.hppでは、BondLengthInteractionなどが、使うポテンシャルの名前を読み、read_potential.hppで定義されている関数を呼び出しています。

ではまずread_interaction.hppで、BondLengthInteractionとQuarticPotentialの組み合わせを読み込めるようにしましょう。 mjolnir/input/read_interaction.hppを好きなエディタで開いてみてください。最初に様々なポテンシャルをincludeしていますが、コード部分に入るとすぐに、以下のような関数が見つかると思います。

template<typename traitsT>
std::unique_ptr<LocalInteractionBase<traitsT>>
read_bond_length_interaction(
    const typename LocalInteractionBase<traitsT>::connection_kind_type kind,
    const toml::Table& local)
{
    // TOMLテーブル内のpotential名を取り出す
    const auto potential = toml::get<std::string>(
            toml_value_at(local, "potential", "[forcefield.local]"));

    if(potential == "Harmonic")
    {
        using potential_t = HarmonicPotential<traitsT>;
        // read_potential.hppで定義されているread_harmonic_potentialを呼び出す
        return make_unique<BondLengthInteraction<traitsT, potential_t>>(
            kind, read_harmonic_potential<traitsT, 2>(local));
    }
    else if
    ...

この関数はTOMLファイルのテーブル（toml::Table）を受け取って、ポテンシャルの名前を受け取り、適切な関数を読んでいます。なので、追加したいポテンシャルの名前を決め（今回はQuarticです）、同様の分岐をelse ifの中に付け加えてください。

続いて、read_potential.hppで実際にパラメータを読む関数を付け足しましょう。

To be written