MDプログラムの作成ガイド

NB-LIBの目的は、研究者が分子シミュレーションをプログラム的に定義できるようにすることです。従来、これらは、複数の実行可能ファイルと、手動で定義されたワークフロー、そして「ブラックボックス」シミュレーションエンジンを使用して行われてきました。NB-LIBは、ユーザーがより詳細なレベルで、さまざまな新しいシミュレーションおよび分析ワークフローをスクリプト化することを可能にします。

NB-LIBが提供する柔軟性により、さまざまなユースケースを実現できます。これには、カスタムの更新ルール、カスタムフォースの定義、シミュレーション群のオーケストレーションなどが含まれます。また、NB-LIBは、従来のMDシミュレーションと分析の作成も可能にします。

このドキュメントでは、NB-LIBのAPIを使用して、|Gromacs|パッケージに含まれる機能を利用してMDプログラムを作成する手順について説明します。

グローバル定義

NB-LIB プログラムは C++ で記述されるため、I/O や高度なタスクに関連するヘッダーファイルを含める必要があります。さらに、NB-LIB が提供するさまざまな機能や抽象化に対応するためのヘッダーファイルも含める必要があります。これらのヘッダーファイルは、こちらから直接コピーできます。最後に、ライブラリで定義されたデータ構造には、nblib という名前空間を使用します。ブロックの最後の行を使用すると、関数やデータ構造が使用されるたびに、この指定を省略できます。

#include <cstdio>

#include "nblib/box.h"
#include "nblib/forcecalculator.h"
#include "nblib/integrator.h"
#include "nblib/molecules.h"
#include "nblib/nbkerneloptions.h"
#include "nblib/particletype.h"
#include "nblib/simulationstate.h"
#include "nblib/topology.h"

using namespace nblib;

粒子のデータを定義する

// Parameters from a GROMOS compatible force-field 2016H66

struct OWaterAtom
{
    ParticleName         name = "Ow";
    Mass                 mass = 15.999;
    C6                   c6   = 0.0026173456;
    C12                  c12  = 2.634129e-06;
};

struct HwAtom
{
    ParticleName         name = "Hw";
    Mass                 mass = 1.00784;
    C6                   c6   = 0.0;
    C12                  c12  = 0.0;
};

struct CMethAtom
{
    ParticleName         name = "Cm";
    Mass                 mass = 12.0107;
    C6                   c6   = 0.01317904;
    C12                  c12  = 34.363044e-06;
};

struct HcAtom
{
    ParticleName         name = "Hc";
    Mass                 mass = 1.00784;
    C6                   c6   = 8.464e-05;
    C12                  c12  = 15.129e-09;
};

この種の構造は、システム内の粒子タイプ数と同じ数だけ存在できます。 データをこのように整理することは厳密には必要ありませんが、明確にするために示されています。 上記のように、単一の要素に対応する複数の粒子が存在することがあり、分子の文脈によって原子質量が異なる場合があります。 例えば、カルボキシル基中の炭素原子は、メチル基中の炭素原子とは異なるパラメータを持つ可能性があります。 標準的な力場からパラメータセットを取得するか、新しい化合物や力場を研究するために新しいパラメータを生成できます。 この例は、2016H66 パラメータセット からのものです。

座標、速度と力のバッファの定義

std::vector<gmx::RVec> coordinates = {
    { 0.794, 1.439, 0.610 }, { 1.397, 0.673, 1.916 }, { 0.659, 1.080, 0.573 },
    { 1.105, 0.090, 3.431 }, { 1.741, 1.291, 3.432 }, { 1.936, 1.441, 5.873 },
    { 0.960, 2.246, 1.659 }, { 0.382, 3.023, 2.793 }, { 0.053, 4.857, 4.242 },
    { 2.655, 5.057, 2.211 }, { 4.114, 0.737, 0.614 }, { 5.977, 5.104, 5.217 },
};

std::vector<gmx::RVec> velocities = {
    { 0.0055, -0.1400, 0.2127 }, { 0.0930, -0.0160, -0.0086 }, { 0.1678, 0.2476, -0.0660 },
    { 0.1591, -0.0934, -0.0835 }, { -0.0317, 0.0573, 0.1453 }, { 0.0597, 0.0013, -0.0462 },
    { 0.0484, -0.0357, 0.0168 }, { 0.0530, 0.0295, -0.2694 }, { -0.0550, -0.0896, 0.0494 },
    { -0.0799, -0.2534, -0.0079 }, { 0.0436, -0.1557, 0.1849 }, { -0.0214, 0.0446, 0.0758},
};

std::vector<gmx::RVec> forces = {
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
};

std::vector を使用して、gmx::RVec の 3D ベクトルデータを格納するための特定のデータ型を利用して、粒子の座標を初期化できます。RVec に関する Doxygen ページはこちら.

MDプログラムの作成

同様に、あらゆる基本的なC++プログラムと同様に、``main()``関数が必要です。

パーティクルタイプの定義

int main()
{
    // Bring the parameter structs to scope
    OwAtom      owAtom;
    HwAtom      hwAtom;
    CMethAtom   cmethAtom;
    HcAtom      hcAtom;

    // Create the particles
    ParticleType Ow(owAtom.name, owAtom.mass);
    ParticleType Hw(hwAtom.name, hwAtom.mass);
    ParticleType Cm(cmethAtom.name, cmethAtom.mass);
    ParticleType Hc(hcAtom.name, hcAtom.mass);

上記と同様に、ParticleType のデータ定義に使用されるヘルパー構造は厳密には必要ありませんが、理解を助けるために示されています。ParticleType CMethAtom(ParticleName("Cm"), Mass(12.0107)); の行で十分です。

非結合相互作用を定義する

ParticleTypeInteractions interactions(CombinationRule::Geometric);

// add non-bonded interactions for the particle types
interactions.add(owAtom.name, owAtom.c6, owAtom.c12);
interactions.add(hwAtom.name, hwAtom.c6, hwAtom.c12);
interactions.add(cmethAtom.name, cmethAtom.c6, cmethAtom.c12);
interactions.add(hcAtom.name, hcAtom.c6, hcAtom.c12);

Lennard-Jones の相互作用については、「ParticleTypeInteractions」オブジェクトを定義します。 ``ParticleType」の各粒子は、``C6」および「C12」パラメータに基づいて、他のすべての粒子と相互作用します。 2つの異なる粒子のこれらのパラメータは、「Geometric」または「LorentzBerthelot」の「CombinationRule」を使用して平均されます。 詳細については、こちらを参照してください: [http://manual.gromacs.org/documentation/2019/reference-manual/functions/nonbonded-interactions.html#the-lennard-jones-interaction]. デフォルトでは、「CombinationRule::Geometric」が選択されます。

各粒子の相互作用パラメータを ParticleTypeInteractions オブジェクトに追加します。これにより、すべての ParticleType のペアに対する相互作用が指定されたテーブルが作成されます。次の行列は、CombinationRule::Geometric を使用して作成されたペアごとの C6 パラメータを記述しています。

#	はい	ハードウェア	Cm	Hc
はい	0.0026	0.0	0.42	4.7e-4
ハードウェア	0.0	0.0	0.0	0.0
Cm	0.42	0.0	0.013	1.05e-3
Hc	4.7e-4	0.0	1.05e-3	8.5e-5

特定の相互作用ペアの場合、ユーザーは指定された CombinationRule を、カスタムパラメータで上書きすることも可能です。次のオーバーロードは、Ow と Cm の粒子タイプ間の CombinationRule から計算されたパラメータを置き換えます。

interactions.add("Ow", "Cm", 0.42, 42e-6);

モジュール化された再利用可能なコードを実現するために、複数の ParticleTypeInteractions オブジェクトを組み合わせることができます。 otherInteractions が定義されている場合、これを行うには interactions.merge(otherInteractions) を使用します。

モジュールを定義する

Molecule water("Water");
Molecule methane("Methane");

water.addParticle(ParticleName("O"), Ow);
water.addParticle(ParticleName("H1"), Hw);
water.addParticle(ParticleName("H2"), Hw);

water.addExclusion("H1", "O");
water.addExclusion("H2", "O");

methane.addParticle(ParticleName("C"), Cm);
methane.addParticle(ParticleName("H1"), Hc);
methane.addParticle(ParticleName("H2"), Hc);
methane.addParticle(ParticleName("H3"), Hc);
methane.addParticle(ParticleName("H4"), Hc);

methane.addExclusion("H1", "C");
methane.addExclusion("H2", "C");
methane.addExclusion("H3", "C");
methane.addExclusion("H4", "C");

それでは、分子を構成する粒子を宣言から始めます。文字列識別子は、分子内の特定の粒子を一意に識別する必要があります。また、各粒子に対して、クーロン相互作用の計算のために部分電荷を定義することも可能です。「water.addParticle(ParticleName("O"), Charge(-0.04), Ow);」

非結合相互作用の計算を必要な場合にのみ行うように、除外の設定を行う。たとえば、2つの粒子が結合を共有している場合、結合のポテンシャルエネルギーが非結合項を無視できるようにします。粒子間の自己除外はデフォルトで有効になっています。この設定と、後でリストされるすべての相互作用で使用するために、addParticle() で指定された一意の識別子を使用します。

リスト形式での相互作用を定義する

分子内で、構成要素間の結合、角度、二面角などの相互作用を定義できます。NB-LIBは、一般的に使用される2、3、4中心の相互作用の実装を提供します。

HarmonicBondType ohHarmonicBond(1, 1);
HarmonicBondType hcHarmonicBond(2, 1);

DefaultAngle hohAngle(Degrees(120), 1);
DefaultAngle hchAngle(Degrees(109.5), 1);

//add harmonic bonds for water
water.addInteraction("O", "H1", ohHarmonicBond);
water.addInteraction("O", "H2", ohHarmonicBond);

// add the angle for water
water.addInteraction("H1", "O", "H2", hohAngle);

// add harmonic bonds for methane
methane.addInteraction("H1", "C", hcHarmonicBond);
methane.addInteraction("H2", "C", hcHarmonicBond);
methane.addInteraction("H3", "C", hcHarmonicBond);
methane.addInteraction("H4", "C", hhcHarmonicBondc);

// add the angles for methane
methane.addInteraction("H1", "C", "H2", hchAngle);
methane.addInteraction("H1", "C", "H3", hchAngle);
methane.addInteraction("H1", "C", "H4", hchAngle);
methane.addInteraction("H2", "C", "H3", hchAngle);
methane.addInteraction("H2", "C", "H4", hchAngle);
methane.addInteraction("H3", "C", "H4", hchAngle);

シミュレーションと非結合計算のためのオプションを定義する

// Define a box for the simulation
Box box(6.05449);

// Define options for the non-bonded kernels
NBKernelOptions options;

立方体の境界ボックスは、1つの引数で定義することも、長さ、幅、高さの3つの引数を個別に指定することもできます。

NBKernelOptions は、ハードウェアコンテキストとシミュレーションに必要な関連計算の両方に対するフラグと設定オプションのセットを含んでいます。以下の表は、設定できる可能性のあるオプションを記述しています。

フラグまたは設定オプション	種類	影響
useGpu	ブール値	非結合計算のためにGPUを使用する
numThreads	インテ・ゲ	使用するCPUスレッドの数
nbnxmSimd	列挙	カーネルのSIMD型（SimdAuto/SimdNo/Simd4XM/ Simd2XMM）
useHalfLJOptimization	ブール値	i-clusterの半分のLJ最適化を有効にする
pairlistCutoff	現実	ペアリストとインタラクションの終了点を指定する
computeVirialAndEn ergy	ブール値	エネルギー計算を有効にする
coulombType	列挙	関数（Pme/Cutoff/ReactionField）
useTabulatedEwaldC orr	ブール値	分析的な方法ではなく、表形式のPMEグリッド補正を使用する
numIterations	インテ・ゲ	各カーネルの反復回数を指定する
cyclesPerPair	ブール値	ペア/サイクルではなく、ペア/サイクルを印刷するように設定する
timestep	現実	時間ステップを指定する

トポロジーとシミュレーションの状態を定義する

当社は TopologyBuilder クラスを使用してシステムトポロジーを構築します。 以前に定義した Molecule オブジェクトと ParticleTypesInteractions を、そのパブリック関数を使用して追加します。 buildTopology() 関数を使用して、定義されたエンティティに基づいて構築された、すべての除外、相互作用マップ、およびリストされた相互作用データを含む実際の Topology オブジェクトを取得します。

TopologyBuilder topologyBuilder;

// add molecules
topologyBuilder.addMolecule(water, 10);
topologyBuilder.addMolecule(methane, 10);

// add non-bonded interaction map
topologyBuilder.addParticleTypesInteractions(interactions);

Topology topology = topologyBuilder.buildTopology();

現在、シミュレーションの状態オブジェクトを使用して、システムを完全に記述するために必要なすべての情報を持っています。これは、トポロジー、箱、および粒子座標と速度を使用して構築されます。このオブジェクトは、分析に使用したり、既知の状態からシミュレーションを開始するために使用できる、システムの「スナップショット」として機能します。

SimulationState simulationState(coordinates, velocities, forces, box, topology);

MD ループの作成

これで、システムと問題を完全に記述したので、MDループを作成するために2つのエンティティが必要です。1つ目は「ForceCalculator」、もう1つは「Integrator」です。NB-LIBには「LeapFrog」というインテグレーターが付属していますが、ユーザーは独自のインテグレーターも作成できます。

// The force calculator contains all the data needed to compute forces
ForceCalculator forceCalculator(simulationState, options);

// Integration requires masses, positions, and forces
LeapFrog integrator(simulationState);

// Allocate a force buffer
gmx::ArrayRef<gmx::RVec> userForces(topology.numParticles());

// MD Loop
int numSteps = 100;

for (i = 0; i < numSteps; i++)
{
  userForces = forceCalculator.compute();

  // The forces are not automatically updated in case the user wants to add their own
  std::copy(userForces.begin(), userForces.end(), begin(simulationState.forces()));

  // Integrate with a time step of 1 fs
  integrator.integrate(1.0);
}

return 0;
} // main