OpenMMをステップバイステップで 〜 Part 3:GUIでシミュレーション条件の設定 〜
OpenMMの使い方を順番にたどる記事のPart 3です。Part 1、Part 2の記事でMD計算のためのタンパク質構造の前処理が完了しました。今回の記事ではシミュレーション条件本体のセットアップの流れを追ってみたいと思います。
利用するツールはOpenMMのGUI、OpenMM Setupです。なお、今回の記事では設定のみで計算そのものは実行しません。*1
1. OpenMM Setupでできること
まず最初にOpenMM Setupでできることの概要です。
前々回の記事でPDBファイルを読み込んで構造の前処理の方法を確認しました。この時はPDBの出力で一度終わりとしましたが、そのまま続けてシミュレーション条件を指定するスクリプト(run_openmm_simulation.py
)の作成と保存、さらに実行(Run Simulation
)までを行うことができます。
画面構成はこんな感じ。
上図のように、左側で選んだ条件が右側のスクリプトに逐次的に反映されて設定ファイルが出来上がっていきます。
前処理と同じく指示に従って条件を選択していくだけ!これなら私でもできそう!
2. MD計算の流れと設定項目のおさらい
シミュレーション条件を設定していく前に、MD計算の流れと設定項目のおさらいをしておきましょう。
計算の流れはざっと以下です。
準備直後の系(⓪)は水分子やイオンを適当に配置しただけなので、最初にエネルギー最小化(①)します。ついで各原子にランダムに速度を与え、徐々に系全体を設定温度まで引き上げていきます(平衡化(②))。最後に平衡状態の系においてアンサンブルを達成するようにトラジェクトリを生成します(シミュレーション本番(③))。*2
設定の必要な項目は大体こんな感じ。*3
ごちゃごちゃしてますが、系の初期状態と平衡状態(のアンサンブル)を定めて、そこに至る計算方法(分子力場、インテグレータ)や、状態の制御方法(サーモスタット、バロスタット)を設定します。あとは希望のシミュレーション時間(ステップ数)計算すれば良さそうです。
また、MD計算はヘビーなので計算環境(GPUをつかえるか?)やログ、結果(トラジェクトリ)の出力管理といった項目も大事な設定となってくるようです。
これからOpenMM Setupの4つのタブで順番に設定を確認していきますが、おおよそ以下の対応になっています。
では順番に見ていきましょう!
3. OpenMM Setupでシミュレーション条件を設定しよう!
3-1. 注意!~ 入力手順とデフォルト値の違い
設定を始めていく前に少し注意点があります。OpenMM Setupでは入力の手順で少し表示内容が変わります。
具体的にはPDBファイルを入力とした時、「① 未処理の構造を入力としてOpenMM SetupでClean up
したあと、続けてシミュレーションの設定に進んだ場合(Part 1 記事でご紹介した手順)」と、「② すでに処理済みのPDBファイルを入力として直接シミュレーションの設定に進んだ場合」です。
上図左のように前者(①: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では以下のようになっています。
Nonbonded Method
やConstraints
といった項目があります。右側のスクリプトでは「# System Configuration 以下」に対応していることがわかります。
左側の各欄にカーソルを合わせると「何を設定する項目か?」ヘルプが表示されます。
今回はTeachOpenCADDのT019 · Molecular dynamics simulationやデフォルト値を参考に以下のように設定してみました。
まず、① 長距離静電相互作用(nonbondedMethod
)はPME(粒子・メッシュ・エバルト法、エバルトの方法 - Wikipedia)です。*4
② カットオフ距離(nonbondedCutoff
)はPME法のうち、「粒子部分(実空間における短距離ポテンシャルの直接和)」に関するパラメータで、「直接計算する距離をどの程度の長さにするか?」を定めるカットオフ値のことだと思います。*5
③ Ewald法の誤差許容範囲(ewaldErrorTolerance
)は、「Ewald法で和を取る際に、切り捨て計算することで生じる力の端数の誤差」に、おおよそ等しくなるように設定します。デフォルトでは0.0005
です。*6
④ 拘束条件(constraints
)では、結合長や結合角に関する拘束をオプションで設定することができます。HBonds
とした場合は、「水素原子を含む全ての結合の結合長」を拘束します。他にはAllBonds
やHAngles
といった選択肢があります。*7
⑤ 拘束条件の誤差許容範囲(constraintTolerance
)は、「拘束条件について誤差をどの程度許容するか?」で、1x10-6としています。*8
なお、スクリプト側にある「⑥ 剛体としての水(ridigdWater
)」は水分子の扱いを決めるパラメータです。OpenMMでは水分子を結合長・結合角ともに拘束された、完全に剛直なもの(rigid)として扱います。False
にするとflexible waterとなります。
以上がSystemタブでした。*9
3-3. Integratorタブ
つづいて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タブでした。
なお、上図のスクリプト部分は変数の設定のみです。計算方法の設定は以下のようになっていました。
積分が LangevinMiddleIntegrator
や圧調整が MonteCarloBarostat
となっていることが確認できますね。
3-4. Simulationタブ
続いてSimulationタブです。シミュレーション時間や計算環境の設定を行います。
Google ColabのGPUを利用するつもりで、以下のように設定してみました。
① シミュレーションの長さはトラジェクトリを生成する時間(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
③ プラットフォームでは計算環境を設定します。選択肢にはReference
、CPU
、CUDA
、OpenCL
があります。Google ColabではGPUを使用できるのでCUDAとしています。*15
④ 精度はプラットフォームでCUDA
あるいはOpenCL
を選択した場合に選べるパラメータで、「single
、mixed
、double
」があります。
精度よりも速度を重視するなら「single
(単精度)」で、逆に精度を重視するなら「double
(倍精度)」です。「mixed
」は力をsingle、積分をdoubleで計算する混合モードです。*16
以上がSimulationタブでした。
3-5. Outputタブ
最後にOutputタブです。結果の出力を設定します。色々な形式で出力できて、項目数が多いので分割して見ていきます。
OpenMMではトラジェクトリをPDB、PDBx/mmCIFあるいはDCDフォーマットで保存することができます。*17GUIのOpenMM Setupでは、このうちDCDフォーマットのみを選択できます。*18
まず、① ボックスにチェックを入れて「トラジェクトリを保存する」を選択し、ファイル名と「どれくらいの間隔(step)ごとに結果を出力するか?」を指定します。
上図では「trajectory.dcd」ファイルに「10000 steps (20 ピコ秒)」毎に状態を書き込むことになります。
また、計算のログファイルを残すこともできます(②)。プロダクションタイムのトラジェクトリ以外のさまざまなデータを保存することができます。
保存する内容はData to Write
の中から選べ、エネルギーや速度、温度などがあります。上記では、「log.txt」ファイルに「1000 steps(2 ピコ秒)」毎に書き込みます。*19
①、②で設定した内容はスクリプトではそれぞれ①DCDReporter( )
、②StateDataReporter( )
に引数として反映されています。
OpenMMではシミュレーションの進行具合をチェックポイントファイルというバイナリファイル(.chk
)で残すことができます(③)。サーバーエラーなどでシミュレーションが途中で止まってしまった時に、チェックポイントファイルを読み込んで計算を再開することができるそうです。
上図では、「checkpoint.chk」ファイルに「10000 steps (20 ピコ秒)」毎に上書き保存しています。スクリプトでは③CheckpointReporter( )
です。*20
また、シミュレーションの状態(位置、速度 etc.)をXMLファイルに書き出すこともできます。このXMLファイルを再度読み込むことでシミュレーションを続けることができます。*21
chkには「同じ計算環境、OpenMMのバージョンでしかシミュレーションを再開できない」という欠点がありますが、XMLにはその欠点がありません。一方で、XMLにはサイズが大きくなるという欠点があります。
上図では、シミュレーションの設定(④)と、最後の状態(⑤)をそれぞれXMLで保存しています。
スクリプトでは、④シミュレーションの設定「system.xml、integrator.xml」はテキストファイルに書き込む形式で書かれており、⑤最後の「final_state.xml」はsaveState( )
というメソッドを使った形式になっています。
なお、⑤はXML
以外に、chk
やPDBx/mmCIF
フォーマットを選択することができます。
以上がOutputタブでした。
3-6. スクリプトを保存
これで4つのタブでの設定がすべて終わりました。最後に作成されたスクリプトを保存しておきましょう。
Save Script
をクリックすると「run_openmm_simulation.py」というファイルが保存されます。
4. おわりに
以上、今回はOpenMM Setupを使ったシミュレーション条件の設定をみてみました。GUIでクリックしていくことで逐次的にスクリプトが生成され、同時に確認できるというのはとても分かりやすくて良いですね!
以前の記事と合わせて、タンパク質構造とシミュレーションスクリプトが揃ったことになります。あとはこれが本当に動くのかどうかが心配です・・・。 早速試したいところですが、記事が長くなってしまったので続きは次回に。
記事の途中、設定の説明を試みていますがOpenMM User Guideを適当に訳しただけなので、意味のわからないところ、間違っているところが多いと思います。ご指摘いただけると嬉しいです。
おしまい。
*1:今回の記事はOpenMM User Guideの3. 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
*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