Developers Guide

English

Developer's Guide

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

新規なlocalポテンシャル

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

結合長に関する相互作用はすでに実装されているので（詳細は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<traitsT>&, const real_type) const noexcept {return;}
    const char* name() const noexcept {return "Quartic";}

  private:
    real_type k_, r_0_;
};

あとは、必要なincludeファイル（<mjolnir/core/System.hpp>と、他にもし必要なら<cmath>などです）と、namespace mjolnir、そしてInclude guardを追加して完成です。 このようなクラスの実例には、例えばmjolnir/potential/local/HarmonicPotential.hppを見てください。

ファイル入力

次は、入力ファイルからこのポテンシャルを指定できるようにしましょう。 LocalInteractionは、以下のような入力形式になっています。

[[forcefields.local]]
interaction = "BondLength"
potential = "Quartic" # 今回作ったポテンシャルの名前を決めます。
topology  = "bond"
parameters = [
# 粒子iとjの間にかかるポテンシャルのr0とkを決めます。
{indices = [0, 1], eq = 3.8, k = 10.0},
{indices = [1, 2], eq = 3.8, k = 10.0},
# ...
]

この形式でポテンシャルを読み込めるようにしましょう。 そのためには、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の中に付け加えてください。

    else if(potential == "Quartic")
    {
        using potential_t = QuarticPotential<traitsT>;
        return make_unique<BondLengthInteraction<traitsT, potential_t>>(
            kind, read_quartic_potential<traitsT, 2>(local));
    }

続いて、read_potential.hppで実際にパラメータを読む関数を付け足しましょう。 LocalInteractionは、「ある粒子間に」「あるパラメータで」定義されるので、粒子のインデックスとパラメータのペアを読み込むことになります。 パラメータは先ほど定義したQuarticPotentialクラスそのものに渡されます。 ここで粒子のインデックスは、結合長ポテンシャルなら2つ、結合角ポテンシャルなら3つ、二面角ポテンシャルなら4つと変化するので、std::array<std::size_t, N>で表現されます。このNはread_interactionで渡されています。

        return make_unique<BondLengthInteraction<traitsT, potential_t>>(
            kind, read_quartic_potential<traitsT, 2>(local));
//                                                ^

よって入力される情報の単位は、std::pair<std::array<std::size_t, N>, QuarticPotential<traitsT>>ということになります。 さらに、LocalInteractionは通常、多くの粒子のペアの間に定義されるものなので、read_quartic_potentialはstd::vector<std::pair<...>>を返します。なので、関数の型は以下のようになります。

template<typename traitsT, std::size_t N>
std::vector<std::pair<std::array<std::size_t, N>, QuarticPotential<traitsT>>>
read_quartic_potential(const toml::Table& local)
{
    ...
}

この返り値が以下の入力に対応していることはもうお気付きでしょう。

parameters = [
# 粒子番号とパラメータ
{indices = [0, 1], eq = 3.8, k = 10.0},
# ...の配列
]

では、中身を書いていきましょう。まず、TOMLファイルのparametersという変数を配列として取り出します。

template<typename traitsT, std::size_t N>
std::vector<std::pair<std::array<std::size_t, N>, QuarticPotential<traitsT>>>
read_quartic_potential(const toml::Table& local)
{
    const auto& params = toml_value_at(local, "parameters",
        "[forcefield.local] for Quartic potential"
        ).cast<toml::value_t::Array>();
}

そして、その中身を読んでstd::vectorに格納し、返します。parametersの一つ一つの要素はTOMLのインラインテーブルなので、それぞれTableに変換し、そこからtoml::getを使って必要な値を取り出してstd::vectorにpush_backしていけば済みます。

template<typename traitsT, std::size_t N>
std::vector<std::pair<std::array<std::size_t, N>, QuarticPotential<traitsT>>>
read_quartic_potential(const toml::Table& local)
{
    const auto& params = toml_value_at(local, "parameters",
        "[forcefield.local] for Quartic potential"
        ).cast<toml::value_t::Array>();

    std::vector<std::pair<
        std::array<std::size_t, N>, QuarticPotential<traitsT>
        >> retval;

    for(const auto& item : params)
    {
        using indices_t = std::array<std::size_t, N>;
        using real_type = typename traitsT::real_type;

        const auto& parameter = item.cast<toml::value_t::Table>();
        auto indices = toml::get<indices_t>(toml_value_at(parameter, "indices",
                "element of [[parameters]] in quartic potential"));
        auto r0 = toml::get<real_type>(toml_value_at(parameter, "eq",
                "element of [[parameters]] in quartic potential"));
        auto k  = toml::get<real_type>(toml_value_at(parameter, "k",
                "element of [[parameters]] in quartic potential"));

        retval.push_back(std::make_pair(indices, QuarticPotential<traitsT>(k, r0)));
    }
    return retval;
}

これで終わりです。新しい結合長ポテンシャルが使えるようになりました。

さらに、ここまでは言及しませんでしたが、read_quartic_potentialは結合長のみならず結合角、二面角相互作用からもそのまま呼び出せます。唯一の違いはかかる粒子の数が異なることですが、それはread_interactionからread_quartic_potentialを呼ぶ時のtemplate引数Nで指定されるので、read_quartic_potentialを変更する必要はありません。追加したい場合は、今回結合長について行ったのと同じく、結合角、二面角相互作用に対してelse if(potential == "Quartic")節を追加するだけで完了です。

なので、タイトルでは「新規な結合長ポテンシャルの追加」と言っていましたが、これは実質新規なポテンシャルをLocal相互作用の全てに追加する方法です。Local相互作用のためのポテンシャルを定義して読み込めるようにするだけで、自然に全てのLocal相互作用から使えるようになっているのが、Mjolnirの特徴です。