magattacaのブログ

日付以外誤報

OpenMMをステップバイステップで 〜 Part 3:GUIでシミュレーション条件の設定 〜

OpenMMの使い方を順番にたどる記事のPart 3です。Part 1Part 2の記事でMD計算のためのタンパク質構造の前処理が完了しました。今回の記事ではシミュレーション条件本体のセットアップの流れを追ってみたいと思います。

利用するツールはOpenMMGUIOpenMM Setupです。なお、今回の記事では設定のみで計算そのものは実行しません。*1

1. OpenMM Setupでできること

まず最初にOpenMM Setupでできることの概要です。

前々回の記事でPDBファイルを読み込んで構造の前処理の方法を確認しました。この時はPDBの出力で一度終わりとしましたが、そのまま続けてシミュレーション条件を指定するスクリプトrun_openmm_simulation.py)の作成保存、さらに実行Run Simulation)までを行うことができます。

画面構成はこんな感じ。

f:id:magattaca:20220402183941p:plain
シミュレーション条件設定の画面構成

上図のように、左側で選んだ条件が右側のスクリプト逐次的に反映されて設定ファイルが出来上がっていきます。

前処理と同じく指示に従って条件を選択していくだけ!これなら私でもできそう!

2. MD計算の流れと設定項目のおさらい

シミュレーション条件を設定していく前に、MD計算の流れと設定項目のおさらいをしておきましょう。

計算の流れはざっと以下です。

f:id:magattaca:20220402184012p:plain
MD計算のおおまかな流れ

準備直後の(⓪)は水分子やイオンを適当に配置しただけなので、最初にエネルギー最小化(①)します。ついで各原子にランダムに速度を与え、徐々に系全体を設定温度まで引き上げていきます(平衡化(②))。最後に平衡状態の系においてアンサンブルを達成するようにトラジェクトリを生成します(シミュレーション本番(③))。*2

設定の必要な項目は大体こんな感じ。*3

f:id:magattaca:20220402184033p:plain
各工程で設定が必要となる項目

ごちゃごちゃしてますが、系の初期状態と平衡状態(のアンサンブル)を定めて、そこに至る計算方法分子力場インテグレータ)や、状態の制御方法サーモスタットバロスタット)を設定します。あとは希望のシミュレーション時間(ステップ数)計算すれば良さそうです。

また、MD計算はヘビーなので計算環境GPUをつかえるか?)やログ、結果(トラジェクトリ)の出力管理といった項目も大事な設定となってくるようです。

これからOpenMM Setup4つのタブで順番に設定を確認していきますが、おおよそ以下の対応になっています。

f:id:magattaca:20220402184052p:plain
各タブと設定内容の対応

では順番に見ていきましょう!

3. OpenMM Setupでシミュレーション条件を設定しよう!

3-1. 注意!~ 入力手順とデフォルト値の違い

設定を始めていく前に少し注意点があります。OpenMM Setupでは入力の手順で少し表示内容が変わります。

具体的にはPDBファイルを入力とした時、「① 未処理の構造を入力としてOpenMM SetupでClean upしたあと、続けてシミュレーションの設定に進んだ場合(Part 1 記事でご紹介した手順)」と、「② すでに処理済みPDBファイルを入力として直接シミュレーションの設定に進んだ場合」です。

f:id:magattaca:20220402184128p:plain
手順によって初期パラメータが異なる

上図左のように前者(①:Yes. Let's clean it up now.)では、シミュレーションの設定項目のほとんどが空欄でひとつずつ自分で選んでいきます。一方、後者(②:No. my file is all ready to simulate.)ではデフォルトの設定(数値)が最初から記入されていて、変更したい箇所だけを変えればOKです。

Clean upをした後でデフォルトの数値設定が知りたい場合は、処理後のPDBファイルを一度保存して、再度読み込んで最初からはじめればOKです。

3-2. Systemタブ

それではひとつずつタブを見ていきましょう。まずはSystemタブです。

「System(系)」とついていますが、ここでは主に長距離相互作用静電相互作用)の計算方法を設定します。

先の項(3-1.)で見た通り、タンパク質構造の前処理からそのまま進んだ時は項目の大半が空欄です。Systemでは以下のようになっています。

f:id:magattaca:20220402184201p:plain
Systemタブ(初期値がほぼ空白の場合)

Nonbonded MethodConstraintsといった項目があります。右側のスクリプトでは「# System Configuration 以下」に対応していることがわかります。

左側の各欄にカーソルを合わせると「何を設定する項目か?」ヘルプが表示されます。

f:id:magattaca:20220402184223p:plain

今回はTeachOpenCADDT019 · Molecular dynamics simulationやデフォルト値を参考に以下のように設定してみました。

f:id:magattaca:20220402184245p:plain
Systemの条件設定例

まず、① 長距離静電相互作用(nonbondedMethod)はPME(粒子・メッシュ・エバルト法、エバルトの方法 - Wikipedia)です。*4

カットオフ距離nonbondedCutoff)はPME法のうち、「粒子部分(実空間における短距離ポテンシャルの直接和)」に関するパラメータで、「直接計算する距離をどの程度の長さにするか?」を定めるカットオフ値のことだと思います。*5

Ewald法の誤差許容範囲ewaldErrorTolerance)は、「Ewald法で和を取る際に、切り捨て計算することで生じる力の端数の誤差」に、おおよそ等しくなるように設定します。デフォルトでは0.0005です。*6

拘束条件constraints)では、結合長や結合角に関する拘束をオプションで設定することができます。HBondsとした場合は、「水素原子を含む全ての結合の結合長」を拘束します。他にはAllBondsHAnglesといった選択肢があります。*7

拘束条件の誤差許容範囲constraintTolerance)は、「拘束条件について誤差をどの程度許容するか?」で、1x10-6としています。*8

なお、スクリプト側にある「⑥ 剛体としての水ridigdWater)」は水分子の扱いを決めるパラメータです。OpenMMでは水分子を結合長・結合角ともに拘束された、完全に剛直なもの(rigid)として扱います。Falseにするとflexible waterとなります。

以上がSystemタブでした。*9

3-3. Integratorタブ

つづいてIntegratorタブです。積分計算やアンサンブルについて設定します。

こんな感じで設定してみました。

f:id:magattaca:20220402184313p:plain
Integratorタブの設定例

まず、① 積分計算の時間間隔Step Size, dt)は「0.002ピコ秒(2フェムト秒)」としました。

統計的アンサンブルStatistical Ensemble)には、「Constant pressure, temperature, NPT」と「Constant volume, temperature, NVT」の選択肢があります。今回はNPTとしました。

尚、GUIのOpenMM Setupでは「アンサンブルを達成するための計算方法(温度や圧力を一定に保つための計算方法)」を選択する項目は無いようです。NPTを選択するとデフォルトで「温度LangevinMiddleIntegrator圧力MonteCarloBarostat」が設定され、スクリプトに反映されます。

温度temperature)は「300K(27℃)」にしました。

摩擦係数Friction Coefficient, friction)は温度制御に関わるパラメータです。デフォルトのLangevinMiddleIntegratorはランジュバン動力学に基づくものなので、ランジュバン方程式の摩擦項の係数をパラメータとして設定します。*10*11

圧力pressure)の単位は「atm(標準大気圧)」です。

圧力調整の間隔Barostat Interval, barostatInterval)は、「何ステップごとに圧力調整を行うか?」です。①で「Step size : 0.002 ps」としているので、⑥で間隔を「25 steps」とすると「0.002 x 25 = 0.050 ps」毎に調整を行います。*12*13

以上がIntegratorタブでした。

なお、上図のスクリプト部分は変数の設定のみです。計算方法の設定は以下のようになっていました。

f:id:magattaca:20220402184339p:plain
積分法、圧調整法のスクリプト

積分LangevinMiddleIntegratorや圧調整が MonteCarloBarostat となっていることが確認できますね。

3-4. Simulationタブ

続いてSimulationタブです。シミュレーション時間や計算環境の設定を行います。

Google ColabのGPUを利用するつもりで、以下のように設定してみました。

f:id:magattaca:20220402184357p:plain
Simulationタブの設定例

シミュレーションの長さはトラジェクトリを生成する時間(production)で、単位はstepです。

Integratorタブで積分Step Sizeを「0.002 ps」としました。従って、上図のように「500,000 steps」とすると、「0.002 x 500,000 = 1,000 ps = 1 ns」となり、1 ナノ秒のシミュレーションとなります。

平衡化の長さも同様にstepで長さを指定します。上図では「50,000 steps」(0.1 ナノ秒)としました。

尚、上記の時間設定はかなり適当なので参考にしないでください。目的にあった適切な値を選択してください。*14

プラットフォームでは計算環境を設定します。選択肢にはReferenceCPUCUDAOpenCLがあります。Google ColabではGPUを使用できるのでCUDAとしています。*15

精度はプラットフォームでCUDAあるいはOpenCLを選択した場合に選べるパラメータで、「singlemixeddouble」があります。

精度よりも速度を重視するなら「single(単精度)」で、逆に精度を重視するなら「double(倍精度)」です。「mixed」はをsingle、積分をdoubleで計算する混合モードです。*16

以上がSimulationタブでした。

3-5. Outputタブ

最後にOutputタブです。結果の出力を設定します。色々な形式で出力できて、項目数が多いので分割して見ていきます。

f:id:magattaca:20220402184549p:plain
Outpuタブの設定例(1)

OpenMMではトラジェクトリをPDBPDBx/mmCIFあるいはDCDフォーマットで保存することができます。*17GUIのOpenMM Setupでは、このうちDCDフォーマットのみを選択できます。*18

まず、① ボックスにチェックを入れて「トラジェクトリを保存する」を選択し、ファイル名と「どれくらいの間隔(step)ごとに結果を出力するか?」を指定します。

上図では「trajectory.dcd」ファイルに「10000 steps (20 ピコ秒)」毎に状態を書き込むことになります。

また、計算のログファイルを残すこともできます(②)。プロダクションタイムのトラジェクトリ以外のさまざまなデータを保存することができます。

保存する内容はData to Writeの中から選べ、エネルギーや速度、温度などがあります。上記では、「log.txt」ファイルに「1000 steps(2 ピコ秒)」毎に書き込みます。*19

①、②で設定した内容はスクリプトではそれぞれ①DCDReporter( )、②StateDataReporter( )に引数として反映されています。

f:id:magattaca:20220402184612p:plain
Outpuタブの設定例(2)

OpenMMではシミュレーションの進行具合をチェックポイントファイルというバイナリファイル(.chk)で残すことができます(③)。サーバーエラーなどでシミュレーションが途中で止まってしまった時に、チェックポイントファイルを読み込んで計算を再開することができるそうです。

上図では、「checkpoint.chk」ファイルに「10000 steps (20 ピコ秒)」毎に上書き保存しています。スクリプトでは③CheckpointReporter( )です。*20

f:id:magattaca:20220402184632p:plain
Outpuタブの設定例(3)

また、シミュレーションの状態(位置、速度 etc.)をXMLファイルに書き出すこともできます。このXMLファイルを再度読み込むことでシミュレーションを続けることができます。*21

chkには「同じ計算環境、OpenMMのバージョンでしかシミュレーションを再開できない」という欠点がありますが、XMLにはその欠点がありません。一方で、XMLにはサイズが大きくなるという欠点があります。

上図では、シミュレーションの設定(④)と、最後の状態(⑤)をそれぞれXMLで保存しています。

スクリプトでは、④シミュレーションの設定「system.xmlintegrator.xml」はテキストファイルに書き込む形式で書かれており、⑤最後の「final_state.xml」はsaveState( )というメソッドを使った形式になっています。

なお、⑤はXML以外に、chkPDBx/mmCIFフォーマットを選択することができます。

以上がOutputタブでした。

3-6. スクリプトを保存

これで4つのタブでの設定がすべて終わりました。最後に作成されたスクリプトを保存しておきましょう。

f:id:magattaca:20220402184658p:plain
スクリプトGet!

Save Scriptをクリックすると「run_openmm_simulation.py」というファイルが保存されます。

4. おわりに

以上、今回はOpenMM Setupを使ったシミュレーション条件の設定をみてみました。GUIでクリックしていくことで逐次的にスクリプトが生成され、同時に確認できるというのはとても分かりやすくて良いですね!

以前の記事と合わせて、タンパク質構造とシミュレーションスクリプトが揃ったことになります。あとはこれが本当に動くのかどうかが心配です・・・。 早速試したいところですが、記事が長くなってしまったので続きは次回に。

記事の途中、設定の説明を試みていますがOpenMM User Guideを適当に訳しただけなので、意味のわからないところ、間違っているところが多いと思います。ご指摘いただけると嬉しいです。

おしまい。

*1:今回の記事はOpenMM User Guide3. Runnning Simulationの記述を参考にしています。

*2:参考:MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜

*3:イラストはいらすとやさんより

*4:プルダウンで選べるPME以外の選択肢は図の右上

*5:OpenMM User Guide 3.1. A First Exampleでは「cutoff for direct space interactions」となっています。湯どうふさんのYouTubeの解説動画「分子シミュレーション実行時のクーロンポテンシャルの効率的な計算法:Ewald法とParticle Mesh Ewald法【理論化学、計算科学】」がとても分かりやすいです。

*6:OpenMM User Guide 3.6.5. Nonbonded Interactions

*7:OpenMM User Guide 3.6.6. Constraints

*8:OpneMMのGitHub Constraint Tolerance Guidelinesでは「1x10-5と1x10-6のどちらが良いか?」ディスカッションされています。「温度一定なら1x10-5、エネルギー一定なら1x10-6」といった意見が出ています。

*9:今回は設定しませんでしたがHydogen Mass Repartitioningは「重原子に結合する水素原子の重さを変えて、動きの速い水素原子を遅くする」ものだそうです。OpenMM User Guide 3.6.7. Heacy Hydrogens

*10:ランジュバン動力学 - wikipedia

*11:OpenMM User Guide 3.6.8. Integrators

*12:OpenMM User Guide 3.6.10. Pressure Coupling

*13:MDの圧力調整におけるモンテカルロ法(Monte Carlo barostat)については、森野キートスさんのブログ記事「AmberでMD計算するときの圧力制御法の違い」などが参考になります。

*14:OpenMM Setupのデフォルト値では、「Step Size :0.004 ps、Simulation Length :1000000 steps、Equilibration Length:1000 steps」のようです。

*15:それぞれの選択肢の詳細はOpenMM User Guide 7.7. Platformsをご参照ください。Referenceはパフォーマンスは考えずに読みやすいコードにしたもので、他の選択肢はそれぞれのプラットフォームに合わせてパフォーマンスを最適化しています。

*16:各プラットフォームで選択できる精度を含む個別のパラメータについてはOpenMM User Guide 10. Platform-Specific Propertiesをご参照ください。

*17:OpenMM User Guide 3.6.13. Writing Trajectories

*18:DCDはMD計算のトラジェクトリでよく使われるフォーマットのようです。参考:DCD Plugin, Version 1.11

*19:OpenMM User Guide 3.6.14. Recording Other Data

*20:OpenMM User Guide 3.6.15. Saving Simulation Progress and Results

*21:OpenMM User Guide 3.6.15. Saving Simulation Progress and Results