-
Notifications
You must be signed in to change notification settings - Fork 7
Design
Mjolnirの開発における目標は、大きく分けて以下のようなものがあります。
- 柔軟性: 新規な力場を簡単に実装できるようにすること。
- 透明性: ユーザーに指定されたことのみを行うこと。
- 高速性: 可能な限り、高速に動作すること。
これらの特徴はどれも、シミュレーション研究を行う上で重要なものです。柔軟で簡潔なコードは実装にかける時間を短縮し、計算機実験と考察に費やせる時間を延ばすことで、科学をより促進するでしょう。ユーザーが指定したことだけを行うことで、内部で何が起きているかが明白になり、想定していなかったアーティフアクトを回避できるでしょう。高速に動くコードはサンプリング効率を上げ、仮説・検証・考察のサイクルを早め、知見の蓄積を助けるでしょう。
これらの目標を同時に達成するため、Mjolnirではいくつかの設計上の方針を持っています。
Mjolnirは、互いに独立した複数の構成要素の組み合わせからできています。構成要素にもそれぞれのカテゴリがあり、同じカテゴリに属するものは同じ使われ方をします。そして同じカテゴリに属するそれぞれの構成要素はお互いの存在を知らず、相互に全く依存しません。
このカテゴリには、例えば「時間積分」、「分子間相互作用」、「ポテンシャル関数」などがあります。Mjolnirのシミュレーションは、それぞれのカテゴリに属するものから適切なものを選択した組み合わせとして実現されます。一部を共有する組み合わせを実装するときに、重複したコードを書く必要は全くありません。
この設計により、新規な分子力場の実装が飛躍的に容易になります。それぞれの構成要素は互いのことを全く知らないので、新たに全く知らない構成要素が追加されても、影響を全く受けません。新たな分子力場を実装するのに必要なのは、新しい構成要素の追加と、それを使った新しい組み合わせを設定ファイルから読み込めるようにすることだけです。
この方針は、一言で言うと「疎結合な設計をできるだけ完璧に実践する」ということでもあります。
Mjolnirは、できるだけtemplateを使って構成要素の組み合わせを実現します。そして構成要素の複数の組み合わせをコンパイル時に生成します。
この方針により、人間が書くべきコードの量が減少します。人間が書いた個々の構成要素だけから、コンパイラが膨大な組み合わせを生成するからです。そして、コンパイラにあり得る組み合わせを教えておくことによって、コンパイラが個別の組み合わせに対してコードを最適化することが可能になります。
また、派生型による仮想関数のオーバーライドを行うのは、基本的に純粋仮想関数のみを対象とします。さらに、派生は1回きりとし、一度派生した型から更に派生はしません。つまり、Mjolnirにおいては継承はインターフェース定義以外の文脈で使いません。
この方針は、一言で言うと「コンパイル時にできることはコンパイル時にする」ということでもあります。
原則としてMjolnirは分離できる部分を分離しますし、コードの重複を避けます。ただし、非常に特殊な場合においては、限定的にこれらの制限を解除します。ただし、原則に従わないコードを書く場合は、多少の議論を要するでしょう。
例えば、非常に特殊な分子力場の中には、ポテンシャル関数を変更するというアイデアが意味をなさないものがあります。そのような場合、ポテンシャル関数を無理に分離はしません。これは、コードに無用な抽象化や複雑性を導入しないようにするための理由からです。「無用な」というのは、「異なる組み合わせで使われないことが明らかな部分を、労力を注いで異なる組み合わせで使えるようにする」状況を意味しています。
他の場合として、最適化があります。環境依存の高度な最適化を行う場合、通常時に使用されるコードに導入しようとすると、複数の環境依存マクロや定数による分岐が追加されるでしょう。これは、それ以外の環境を想定している開発者にとってはノイズになり、ロジックへの集中を妨げるでしょう。そのような場合、環境依存のコードを全く別に作成し、それをまた一つ別の構成要素として取り扱った方が、それぞれのコードがそれぞれの目的に集中できます。
上記のような場合、疎結合で冗長性のない設計から少し離れたところに、明快さや利便性があります。
Mjolnirの目標は明快さと利便性です。ほとんどの場合、疎結合で簡潔なコードがその目標を体現するでしょう。なので、ほぼ常に上述の方針を適用します。しかし、それが機能しない例外的なケースでは、本来の目標を優先します。
MjolnirはScienceのためのパッケージですが、Software Engineeringを一切軽視しません。 Software Engineeringとは、向き合っている問題を深く知り、それを適切なレベルに抽象化して、過不足なく表現することです。それに本気で取り組み実践していくことが、結果的にScienceを加速していくと、私は信じて開発を進めています。
以下では代表的なクラスの設計を説明します。コードを見る前に読んでおくと、理解が速くなるでしょう。
系はSystemというクラスによって表現されます。このクラスは、全てのParticleに関する情報、つまり質量、位置、速度、力、そして系全体に関わる情報（例えばアンサンブルが緩和するべき温度・圧力や、境界条件など）を格納しているクラスです。また、後述しますが、ForceFieldから生成されたTopologyクラスも格納しています。
このクラスの責任は、その系の情報をあまさず格納し、必要に応じて引き出せるようにすることです。
Systemは、境界条件による細かな差異を吸収する働きも持っています。Systemに2点間の距離を問い合わせれば、境界条件に応じた適切な距離を返すようになっています。例えば、周期境界条件を採用しているなら、境界のそばにいる粒子間の最短距離は境界をまたぐ経路かもしれません。この関数はそれを確認し、最小距離やその経路に対応するベクトルを計算します。
そのような関数は、例えば境界が存在せず系が無限に広がっているとき、するべきことがありません。そのような場合にコンパイラが関数呼び出しを削除する機会を設けるため、境界条件は関数ポインタや継承ではなくtemplateを使って実装されています。そのようにすることで、実行時ではなくコンパイル時に呼び出される関数が決まり、コンパイラが何もしない関数を呼び出すオーバーヘッドを削除することが可能になります。
BoundaryConditionはクラスとして提供されます。このクラスは、境界の形状と位置を格納し、そしてそれが粒子に与える影響を計算することです。
このクラスの責任は、境界の形状と位置を格納し、粒子間の距離や粒子の位置を適切なものに変更することです。
分子力場は、ForceField クラスによって表現されます。このクラスは、どの粒子間に、どのような関数形の、どのような相互作用がかかるかを格納しています。例えば、「粒子0と1の間に」、「調和振動子のようなポテンシャル関数の」、「結合距離に関する相互作用がかかる」、や、「粒子0, 1, 2..., Nの全ての間に」、「距離の逆2乗型の」、「空間中の距離に応じた相互作用がかかる」、というような情報を格納しています。
このクラスの責任は、力場のパラメータを格納し、また粒子にかかる力とエネルギーを計算する方法を提供することです。
ForceFieldには3種類あります。Local, Global, Externalの3つです。
-
LocalForceField- 特定の粒子間にのみ働く力場です。結合距離に応じた調和振動子ポテンシャルや、2面角に応じてかかる
cos型ポテンシャルなどです。
- 特定の粒子間にのみ働く力場です。結合距離に応じた調和振動子ポテンシャルや、2面角に応じてかかる
-
GlobalForceField- 対応するパラメータを持つ粒子間に働く力場です。静電相互作用や分子間力などです。電荷などに対応するパラメータが0、または定義されていない粒子との間には力はかかりません。
-
ExternalForceField- 対応するパラメータを持つ粒子に働く力場です。粒子間ではなく、粒子に働きます。分子を箱の内部に閉じ込めるためのポテンシャルや、特定の分子の位置をある座標に固定する力などです。この力場のみが、系全体のトータルでの並進・回転に影響を与えます。
全てのForceFieldクラスは、実際にはInteractionクラスの集合です。
Interactionクラスは、分子力場の持つ相互作用による力、エネルギーを計算する方法を体現したクラスです。例えば、結合角の値に応じた相互作用が定義されているとき、粒子にかかる力はその粒子の位置の微小変化によってエネルギーがどの程度変化するかによって計算されます。このとき、力の大きさはポテンシャルの形によって決まりますが、力の向きは「結合角に関する相互作用」という部分で決まります。
Interactionクラスの責任は、この力の向きを計算することと、結合角（あるいはそれに相当するパラメータ）を計算することです。
これは何を意味しているかというと、一度結合角にかかるポテンシャルを実装すれば、そのInteractionクラスは他のどのような形の関数にも、例えば調和振動子や、サイン関数や、シグモイドや、ありとあらゆるポテンシャル関数に対してそのまま適用できるということです。
数式で説明したほうがわかりやすいかもしれません。ポテンシャル関数 f が角度 θ に依存するとき、そのポテンシャルから粒子（位置 r ）が受ける力は F = - ∂f(θ)/∂r になります。これは、f の θ による微分と、 θ の位置による微分によって、 F = - ∂f/∂θ・∂θ/∂r のようにわけられます。θ の位置による微分を計算することが、Interactionクラスの仕事の1つです。もう1つの仕事は、エネルギーを計算するときに θ を計算することです。
そして、ポテンシャル関数（と、その微分）を表現しているのがPotentialクラスです。
分子力場の持つ関数は、Potentialクラスによって表されます。このクラスは、ポテンシャルが必要とするパラメータを記憶し、ポテンシャル関数の値、そして微分を計算する方法を提供します。
例えば、LocalInteractionが使うPotentialクラスは、例えばその結合角ポテンシャルに関するパラメータと、ある角度での関数の値、微分の値を計算するメソッドを提供します。LocalInteractionクラスは、このPotentialクラスを必要なだけ持ち、順に使用していきます。
他に、GlobalInteractionが使うPotentialクラスは、そのポテンシャルに関わる粒子が持つパラメータ（電荷など）を格納し、またInteractionクラスが提供するパラメータ（距離など）から、関数の値、微分の値を計算するメソッドを提供します。
ExternalInteractionが使うPotentialクラスも、GlobalInteractionのそれと同様のものです。
以上の構成は、広く使われている相互作用、例えば結合長や結合角、二面角などを複数の異なるポテンシャル関数で使い回すためのものです。このようにすることで、新規なポテンシャル関数を試すときに書かなければならないコードの量は激減します。が、少しばかりデメリットもあります。この節では、そのデメリットを回避するために取った方針について説明します。
まず、多くのポテンシャル関数は非常にシンプルで（調和振動子の場合、k*(x-x0)^2です）、関数呼び出しをするほどでもありません。なので、関数ポインタまたは継承を使って動的にディスパッチせずに、InteractionとPotentialの組み合わせでコンパイル時にコードを生成することで、インライン化の機会を与えることにしました。そのために、PotentialクラスはInteractionクラスのtemplate引数として実装されています。
この方針は、今から説明するもう1つの問題の解決策にもなっています。あるInteractionとPotentialの組み合わせの場合にだけ可能になる最適化がいくつか存在します。この2つを分離する設計では、そのような最適化をする機会はなくなってしまいます。しかし、templateを用いれば、特殊化したクラスを定義することでこの最適化を行うことが可能になります。
もちろん、templateの特殊化は、汎用バージョンのコードと特殊なバージョンのコードがどちらが呼ばれているかの判断を難しくします。そのため、そのような最適化を本当にするかどうかは熟慮（と、各所に十分な量のコメント）を必要とするでしょう。しかし、特殊な最適化を行う選択肢が残っていることは、実行速度の観点からは喜ばしいことです。
唯一残ったデメリットは、エラーメッセージです。コンパイラは複数の組み合わせのコードを生成するので、コードの中のたった一箇所のミスをあり得る組み合わせの数だけ指摘します。よって、1つのミスが数十倍の個数のコンパイルエラーになって現れてきます。これはこの方針上避けようがないデメリットですが、コンパイラがコードを生成するという仕事を正しく果たしている証として前向きに受け取ってもらえるとありがたいです。
場合によっては、非常に特殊な力場があり、「ポテンシャル関数を別のものにする」というアイデアが意味をなさない場合があるでしょう。
そのような場合は、無理にPotentialクラスを分けずに、特製のInteractionを作ることができます。
特定のポテンシャル関数とペアになっていないと意味をなさないInteractionは、単体で全てを実装しましょう。Mjolnirではそれが可能です。
Simulatorクラスは、上記の全てを適切に管理するためのクラスです。これは、ForceFieldやSystem、それとIntegratorを持ち、それらを適切なタイミングで呼び出します。
ForceFieldを途中で切り替えたり、シミュレーション全体を複数並行して実行するようなシミュレーション技法が存在しますが、そのような場合Simulatorクラスは必要な数のForceFieldやSystemを持ち、適切なタイミングで呼び出すことで管理することになります。
このクラスの責任は、他の部分を協調させ、シミュレーションのアルゴリズムを適切に実行することです。
Integratorは、時間発展方法を表現するクラスです。速度Verlet法やLangevin動力学などの数値計算アルゴリズムが対応します。
このクラスの責任は、渡されたSystemとForceFieldを使って、Systemの中の粒子を次の時間ステップまで動かすことです。
以上全てのクラスには、templateパラメータとしてSimulatorTraitsが渡されています。このクラスは、全てのクラスに共通して渡されるtemplateパラメータを一纏めにしたものです。そのようなtemplateパラメータとしては、浮動小数点数の型（float、double等）やベクトル型、境界条件型（BoundaryConditionクラスです）などがあります。また、並列化の方法などもこのクラスによってディスパッチされます。