Developers Guide

Developer's Guide

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

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

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

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

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

ポテンシャルの実装

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

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

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

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

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

    QarticPotential(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 QarticPotential
{
  public:
    typedef typename traitsT::real_type real_type;

    QarticPotential(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 QarticPotential
{
  public:
    typedef typename traitsT::real_type real_type;

    QarticPotential(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 "Qartic";}

  private:
    real_type k_, r_0_;
};

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

ファイル入力

WIP