この記事ではMACEによる簡易的な動力学計算を扱います。この記事を理解すればあなたのPCから手軽に簡易的な分子動力学計算が行えるようになります！ぜひトライしてみてください！
macOS Ventura(13.2.1), Windows 11 version 23H2, Google Colaboratory
自宅でできるin silico創薬の技術書を販売中
新薬探索を試したい方必読！
ITエンジニアである著者の視点から、wetな研究者からもdryの創薬研究をわかりやすく身近に感じられるように解説しています
宣伝
こちらの記事は合成生物学大会iGEMの強豪校であるiGEM-Wasedaさん協力のもと執筆されました。ご協力誠にありがとうございます！
【iGEM-Waseda】は合成生物学の研究を行う早稲田大学の学術サークルです。iGEMと呼ばれる合成生物学の世界大会の世界大会に出場するために日々研究に励んでいらっしゃいます。
本記事では、iGEM2024で日本Undergrad部門で史上初のTOP10に選ばれたプロジェクトの一環として、特にIn Silicoシミュレーションに関わる部分のツールの一部を紹介しています。プロジェクトの詳細については、iGEM-Wasedaの成果報告サイトをご覧いただければ幸いです。
MACEとは
MACEは、機械学習ポテンシャル（Machine Learning Potential）の一種として開発されたツールで、材料内の原子間相互作用を高精度かつ高速に予測できるのが特徴です。従来の量子力学の基本法則に基づいて、物質の性質を予測するための計算手法の第一原理計算（ab initio計算）や分子動力学シミュレーションでは、多大な計算コストを要して原子同士の相互作用エネルギーを評価していました。しかしMACEは、機械学習の手法を用いてより少ない計算コストで高い精度を実現します。
機械学習ポテンシャルとは
機械学習ポテンシャル（MLP）は、コンピュータを使って物質の性質を予測する方法のひとつです。物質の中の原子同士がどのように引き合ったり反発したりするかを計算し、エネルギーや安定性を予測します。従来の第一原理計算という複雑で時間がかかる方法が必要でしたが、MLPでは、機械学習を使ってその計算を効率化しています。たとえば、たくさんの計算データを基に学習することで、原子の動きやエネルギーの変化を高速かつ正確に予測できるようになります。これにより、新しい材料の開発や、既存の材料の改善がスムーズに行えるため、製造や化学、医療分野などで注目されています。
分子動力学計算(MD)とは
分子動力学計算（MD）は、分子や原子の動きをシミュレーションする計算手法で、物理や化学、生物学分野で広く活用されています。MDでは、ニュートンの運動方程式に基づいて各原子の位置と速度を計算し、時間発展を通じて系の挙動を追跡します。主な目的は、微視的なレベルでの分子間相互作用や物質の構造変化、反応経路の理解にあります。計算には多くの原子数が必要とされるため、計算コストが高く、スーパーコンピュータや並列計算技術が用いられることが一般的です。MDは特に、タンパク質の構造変化の予測や、新素材の分子構造の解析、薬剤設計などで重要な役割を果たしています。
MACEによる簡易的な分子動力学計算
早速MACEを使った簡易的な分子動力学計算を行っていきましょう！
ここではMACEを使って、単純な溶媒システム自体がMDシミュレーションで正常に動作していることを評価します。今回は溶媒システムですが、タンパク質を加えた系だと、評価後に RMSDなどの解析を行うことによって、タンパク質の安定性を評価することができます。
この記事では、こちらのページにあるMolecular Dynamics with MACEを改変したコードを実行し、分子動力学計算が実行します。
まず、上記のファイルをGoogle Colaboratoryで開きましょう。
次に、Google ColaboratoryでGPUを使えるようにしておきましょう。ランタイム→ランタイムのタイプを変更→T4 GPUで設定できます。
全コードはこちら
#各種ライブラリのインポート
!pip install mace-torch
from ase.io import read, write
from ase import units
from mace.calculators import mace_off, mace_mp
from ase.md.langevin import Langevin
from ase.md.verlet import VelocityVerlet
from ase.md.andersen import Andersen
from ase.md.velocitydistribution import Stationary, ZeroRotation, MaxwellBoltzmannDistribution
from ase.visualize import view
from ase.optimize import BFGS
import random
import os
import time
import numpy as np
import pylab as pl
from IPython import display
#simpleMD関数定義
def simpleMD(init_conf, T_init, calc, fname, s, T):
"""
init_conf: 初期構造
temp: シミュレーション温度
calc: 計算方法
fname: 出力ファイル名
s: 出力間隔
T: ステップ数
unit: 時間の単位
"""
init_conf.calc = calc
random.seed(701) # ランダムシードの設定
MaxwellBoltzmannDistribution(init_conf, temperature_K=T_init) #初期温度を300Kに設定し、マクスウェル・ボルツマン分布に従って速度を分配
Stationary(init_conf) #系全体の並進運動を除去
ZeroRotation(init_conf) #系全体の回転運動を除去
dyn = Andersen(atoms=init_conf, timestep=1.0*units.fs, temperature_K=T_init, andersen_prob=1e-2)
#ランジュバン動力学を設定。摩擦係数を0.1、温度をtemp、タイムステップをunitに設定。
%matplotlib inline
time_fs = []
temperature = []
energies = []
#remove previously stored trajectory with the same name
os.system('rm -rfv '+fname)
fig, ax = pl.subplots(2, 1, figsize=(6,6), sharex='all', gridspec_kw={'hspace': 0, 'wspace': 0})
def write_frame():
print("Time: {0:.2f} fs, Temperature: {1:.2f} K, Energy: {2:.2f} eV".format(dyn.get_time()/units.fs, dyn.atoms.get_temperature(), dyn.atoms.get_potential_energy()/len(dyn.atoms)))
dyn.atoms.write(fname, append=True)
time_fs.append(dyn.get_time()/units.fs)
temperature.append(dyn.atoms.get_temperature())
energies.append(dyn.atoms.get_potential_energy()/len(dyn.atoms))
ax[0].plot(np.array(time_fs), np.array(energies), color="b")
ax[0].set_ylabel('E (eV/atom)')
# plot the temperature of the system as subplots
ax[1].plot(np.array(time_fs), temperature, color="r")
ax[1].set_ylabel('T (K)')
ax[1].set_xlabel('Time (fs)')
display.clear_output(wait=True)
display.display(pl.gcf())
time.sleep(0.01)
dyn.attach(write_frame, interval=s)
t0 = time.time()
dyn.run(T)
t1 = time.time()
print("MD finished in {0:.2f} minutes!".format((t1-t0)/60))
#ファイルパス、モデルの設定
pdb_path = "solvent_molecs.xyz"
molecule = read(pdb_path)
print(f"molecule: {molecule}")
print(f"molecule.get_positions(): {molecule.get_positions()}")
print(f"molecule.get_chemical_symbols(): {molecule.get_chemical_symbols()}")
mace_calc = mace_off(model='model_swa.model', device="cuda")
#MDの実行
simpleMD(molecule, T_init=1200, calc=mace_calc, fname='mace_md.xyz', s=10, T=2000)各種ライブラリのインポート
ではコードを1行ずつ解説していきます。
以下のコードを実行してMACEに必要なパッケージをインストールします。
!pip install mace-torch以下に示すモジュールもインポートします。
from ase.io import read, write
from ase import units
from mace.calculators import mace_off, mace_mp
from ase.md.langevin import Langevin
from ase.md.verlet import VelocityVerlet
from ase.md.andersen import Andersen
from ase.md.velocitydistribution import Stationary, ZeroRotation, MaxwellBoltzmannDistribution
from ase.visualize import view
from ase.optimize import BFGS
import random
import os
import time
import numpy as np
import pylab as pl
from IPython import display各ライブラリやモジュールの概要は以下の通りです。
ase材料科学や分子シミュレーションに特化したPythonライブラリで、原子構造の操作、分子動力学、構造最適化、計算エンジンとの連携などが可能になります。mace.calculatorsACEライブラリの中で計算器（Calculator）を提供するモジュールです。これを利用することで、MACEの機械学習モデルを用いて原子構造のエネルギーや力を計算し、シミュレーションに組み込むことができます。random乱数を生成するためのライブラリで、シミュレーションの初期条件設定やランダム化が必要な場合に使用します。osファイルやディレクトリの操作、システムコマンドの実行を行うライブラリで、シミュレーション結果のファイル管理などに利用します。time実行時間の計測や遅延の導入に使用するライブラリです。numpy高速な数値計算を可能にするライブラリで、配列や行列の操作、データ解析に便利です。pylabPythonでデータの解析や可視化を行うためのライブラリです。特に、科学計算を目的とした高速な数値計算とプロットを統合的に行う環境を提供します。IPython.displayJupyter Notebook内でリアルタイムにグラフやデータを更新・表示するためのライブラリです。
simpleMD関数定義
次は関数を定義します。関数名はsimpleMDです。返り値なしの関数です。
def simpleMD(init_conf, T_init, calc, fname, s, T):
"""
init_conf: 初期構造
temp: シミュレーション温度
calc: 計算方法
fname: 出力ファイル名
s: 出力間隔
T: ステップ数
unit: 時間の単位
"""
init_conf.calc = calc
random.seed(701) # ランダムシードの設定
MaxwellBoltzmannDistribution(init_conf, temperature_K=T_init) #初期温度を300Kに設定し、マクスウェル・ボルツマン分布に従って速度を分配
Stationary(init_conf) #系全体の並進運動を除去
ZeroRotation(init_conf) #系全体の回転運動を除去
dyn = Andersen(atoms=init_conf, timestep=1.0*units.fs, temperature_K=T_init, andersen_prob=1e-2)
#ランジュバン動力学を設定。摩擦係数を0.1、温度をtemp、タイムステップをunitに設定。
%matplotlib inline
time_fs = []
temperature = []
energies = []
#remove previously stored trajectory with the same name
os.system('rm -rfv '+fname)
fig, ax = pl.subplots(2, 1, figsize=(6,6), sharex='all', gridspec_kw={'hspace': 0, 'wspace': 0})
def write_frame():
print("Time: {0:.2f} fs, Temperature: {1:.2f} K, Energy: {2:.2f} eV".format(dyn.get_time()/units.fs, dyn.atoms.get_temperature(), dyn.atoms.get_potential_energy()/len(dyn.atoms)))
dyn.atoms.write(fname, append=True)
time_fs.append(dyn.get_time()/units.fs)
temperature.append(dyn.atoms.get_temperature())
energies.append(dyn.atoms.get_potential_energy()/len(dyn.atoms))
ax[0].plot(np.array(time_fs), np.array(energies), color="b")
ax[0].set_ylabel('E (eV/atom)')
# plot the temperature of the system as subplots
ax[1].plot(np.array(time_fs), temperature, color="r")
ax[1].set_ylabel('T (K)')
ax[1].set_xlabel('Time (fs)')
display.clear_output(wait=True)
display.display(pl.gcf())
time.sleep(0.01)
dyn.attach(write_frame, interval=s)
t0 = time.time()
dyn.run(T)
t1 = time.time()
print("MD finished in {0:.2f} minutes!".format((t1-t0)/60))simpleMD関数は実際に分子動力学計算(MD)を実行する関数です。引数は以下の通りです。
init_conf: 初期構造を表すAtomsオブジェクトT_init: 初期温度（ケルビン単位）calc: 計算方法（力場または量子化学計算など）fname: シミュレーション結果を保存するファイル名s: 出力間隔（ステップ数）T: シミュレーションの総ステップ数
次に関数内のコードについて解説していきます。
init_conf.calc = calc
random.seed(701) # ランダムシードの設定
MaxwellBoltzmannDistribution(init_conf, temperature_K=T_init) #初期温度を300Kに設定し、マクスウェル・ボルツマン分布に従って速度を分配
Stationary(init_conf) #系全体の並進運動を除去ZeroRotation(init_conf) #系全体の回転運動を除去
ここでは分子動力学計算を行う前の準備を行っています。まず先ランダムシードを固定して再現性を確保しています。次に初期温度を300Kに設定しマクスウェル・ボルツマン分布という確率分布に従って各原子に速度を分配しています。最後に原子全体の並進運動と回転運動を除去します。これで前処理は終了です。
前処理が終わったら次は以下のコードになります。
dyn = Andersen(atoms=init_conf, timestep=1.0*units.fs, temperature_K=T_init, andersen_prob=1e-2)
これはAndersen サーモスタットを実行するコードです。Andersen サーモスタットは、分子動力学 (MD) シミュレーションで温度制御を行う手法の一つです。特に、系を特定の温度に維持しながら、統計力学的に正しい振る舞いを再現するために用いられます。
次はグラフを描画するコードについての解説です。
%matplotlib inline
time_fs = []
temperature = []
energies = []
#remove previously stored trajectory with the same name
os.system('rm -rfv '+fname)
fig, ax = pl.subplots(2, 1, figsize=(6,6), sharex='all', gridspec_kw={'hspace': 0, 'wspace': 0})
def write_frame():
print("Time: {0:.2f} fs, Temperature: {1:.2f} K, Energy: {2:.2f} eV".format(dyn.get_time()/units.fs, dyn.atoms.get_temperature(), dyn.atoms.get_potential_energy()/len(dyn.atoms)))
dyn.atoms.write(fname, append=True)
time_fs.append(dyn.get_time()/units.fs)
temperature.append(dyn.atoms.get_temperature())
energies.append(dyn.atoms.get_potential_energy()/len(dyn.atoms))
ax[0].plot(np.array(time_fs), np.array(energies), color="b")
ax[0].set_ylabel('E (eV/atom)')
# plot the temperature of the system as subplots
ax[1].plot(np.array(time_fs), temperature, color="r")
ax[1].set_ylabel('T (K)')
ax[1].set_xlabel('Time (fs)')
display.clear_output(wait=True)
display.display(pl.gcf())
time.sleep(0.01)%matplotlib inline は、Jupyter Notebook や IPython 環境で使用される マジックコマンド の一つです。このコマンドを実行すると、グラフやプロットがインライン（ノートブックセル内）に直接描画されるようになります。対照的に、このコマンドを使わない場合、グラフは別ウィンドウで表示されることがあります（特にデスクトップ環境で実行した場合）。
その後に3つの配列が定義されています。それぞれの役割は以下の通りです。
time_fs: シミュレーションの時間（フェムト秒）。temperature: 系全体の温度（ケルビン）。energies: 系のポテンシャルエネルギー（eV/原子単位）。
os.system('rm -rfv '+fname)は以前のシミュレーション結果ファイルを削除して、新しいシミュレーション結果を保存する準備をしています。
fig, ax = pl.subplots(2, 1, figsize=(6,6), sharex='all', gridspec_kw={'hspace': 0, 'wspace': 0})このコードは2つのグラフを作成するためのものです。2つのプロットを縦に並べたグラフが作成されます。
figsize=(6, 6): グラフのサイズ（横 6 インチ、縦 6 インチ）。sharex='all': x 軸を共有。gridspec_kw: グラフ間の余白を設定。hspace=0: 縦方向の余白をゼロ。wspace=0: 横方向の余白をゼロ。
次に、関数内関数としてwrite_frameが定義されています。write_frame関数の中身を見ていきましょう。
print("Time: {0:.2f} fs, Temperature: {1:.2f} K, Energy: {2:.2f} eV".format(dyn.get_time()/units.fs, dyn.atoms.get_temperature(), dyn.atoms.get_potential_energy()/len(dyn.atoms)))これは現在のシミュレーション時間、温度、エネルギーを表示するものです。
dyn.get_time()/units.fs: シミュレーション時間をフェムト秒単位で取得。dyn.atoms.get_temperature(): 温度を取得。dyn.atoms.get_potential_energy()/len(dyn.atoms): 系の1原子あたりのポテンシャルエネルギーを計算。
dyn.atoms.write(fname, append=True)これはシミュレーション結果をファイルに保存する処理を行っています。
time_fs.append(dyn.get_time()/units.fs)
temperature.append(dyn.atoms.get_temperature())
energies.append(dyn.atoms.get_potential_energy()/len(dyn.atoms))このコードは取得した時間、温度、エネルギーをそれぞれのリストに追加するためのものです。append 関数は、Python のリストに新しい要素を追加するために使用される関数です。
ax[0].plot(np.array(time_fs), np.array(energies), color="b")
ax[0].set_ylabel('E (eV/atom)')
ax[1].plot(np.array(time_fs), temperature, color="r")
ax[1].set_ylabel('T (K)')
ax[1].set_xlabel('Time (fs)')これは実際に計算した温度とエネルギーをグラフにプロットし、軸ラベルを表示するコードです。
- 上のプロット（
ax[0]）: 時間に対するエネルギー（青）。 - 下のプロット（
ax[1]）: 時間に対する温度（赤）。
display.clear_output(wait=True)
display.display(pl.gcf())
time.sleep(0.01)これはグラフの描画をリアルタイムで行うためのコードです。通常はまず計算をし結果をリスト等に格納してからグラフにプロットするのですが、このコードにより計算しながら同時にプロットをすることができます！
display.clear_output(): 前の表示をクリア。display.display(): 新しいプロットを表示。time.sleep(0.01): 更新間隔を 0.01 秒に設定。
dyn.attach(write_frame, interval=s)
t0 = time.time()
dyn.run(T)
t1 = time.time()
print("MD finished in {0:.2f} minutes!".format((t1-t0)/60))これらのコードは実際の分子動力学計算実行に関係するコードです。attach関数でMD シミュレーション中に結果を記録する関数を指定します。ここでは先ほどのwrite_frame関数を指定しています。intervalはグラフ更新の間隔です。そしてrun関数で実際に分子動力学計算を実行します。引数のTは分子動力学計算を実行する時間を表します。t0とt1は開始時と終了時の時刻を取得して分子動力学計算にかかった時間を計算するために定義されています。
ファイルパス、モデルの設定
関数の定義が終わったらいよいよ分子動力学計算を行う対象ファイルを指定を行います。
まず対象ファイルをダウンロードしましょう。
こちらのページからmodel_swa.modelをダウンロードし、こちらのページからsolvent_molecs.xyzをダウンロードし、プログラムファイルと同一階層のディレクトリに入れてください。solvent_molecs.xyz は溶媒分子の構造や座標を記述したデータファイルです。
この時点で、Google Colaboratory上のフォルダ内は以下のようになっていると思います。
続いて、以下を実行してください。
pdb_path = "solvent_molecs.xyz"
molecule = read(pdb_path)
print(f"molecule: {molecule}")
print(f"molecule.get_positions(): {molecule.get_positions()}")
print(f"molecule.get_chemical_symbols(): {molecule.get_chemical_symbols()}")
mace_calc = mace_off(model='model_swa.model', device="cuda")
pdb_pathに分子動力学計算を行う対象のタンパク質構造ファイルのファイル名を記述してください。タンパク質構造ファイルは様々なファイル形式に対応しています（例: XYZ, CIF, PDB, VASP POSCAR）。この記事ではサンプルデータにあるsolvent_molecs.xyzで試しています。
そしてmace_calcで分子動力学計算の設定を行います。model="model_swa.model”で事前学習モデルの選択を行っています。ここではmodel_swa.modelモデルを選択しています。
最後に、device="cuda”で計算をGPUで行うように設定しています。初めにランタイムのタイプを変更でT4 GPUを選択していれば大丈夫です。このコードを実行すると以下のようなものが表示されると思います。
molecule: Atoms(symbols='CO3C3H2CH8COCOCH2OH2CO3C2H6', pbc=False, molID=...)
molecule.get_positions(): [[-1.10453892 -4.18069029 1.83713901]
[-0.53202987 -3.86270618 2.77864814]
[-0.73509789 -4.11304951 0.57840604]
[-2.13083696 -5.06843948 1.89835703]
[ 0.13405514 -3.05180025 0.36176103]
[-2.91671395 -5.15526438 3.14581203]
[ 0.43207514 -3.09556818 -1.07133293]
[-0.33282685 -2.08445334 0.55586803]
[ 1.1605351 -3.26725221 0.91560405]
[-4.17609406 -5.68575048 2.52319908]
[-2.51292491 -5.5566535 4.09429598]
[-3.09405875 -4.1374445 3.49866199]
[-4.79772091 -6.13905811 3.39235902]
[-4.88089609 -4.98939514 2.20023108]
[-3.95055079 -6.51654911 1.80283999]
[ 1.31437516 -2.41811919 -1.17466593]
[ 0.65923214 -4.11563444 -1.50241697]
[-0.40137285 -2.36518836 -1.57631993]
[-3.00461197 3.96998763 -2.02547503]
[-2.0105288 5.19786978 -2.24847293]
[-2.28143287 6.1203537 -1.41568494]
[-3.19318676 5.59348774 -0.44043297]
[-3.36470985 4.23982096 -0.54081994]
[-3.75257277 4.06419754 -2.75102592]
[-2.58701587 2.96247864 -2.22737694]
[-2.11613488 7.35936356 -1.39925003]
[-4.44593668 4.15878677 -0.46381998]
[-2.87879086 3.74973965 0.31456503]
[ 4.21964693 0.24122973 -0.68176395]
[ 3.34673715 0.74083972 0.17128803]
[ 4.71122122 0.7166627 -1.78589094]
[ 5.02922392 -0.66170031 -0.29166198]
[ 4.2078619 1.82501674 -2.37506008]
[ 4.97652912 -1.0726763 0.964589 ]
[ 3.31383419 1.4959147 -2.97123194]
[ 4.13560104 2.69454861 -1.74539602]
[ 4.91796112 2.03662682 -3.15533805]
[ 5.85860014 -1.77208734 1.13940108]
[ 4.89973307 -0.24120829 1.68790603]
[ 4.03826714 -1.63401127 1.09342408]]
molecule.get_chemical_symbols(): ['C', 'O', 'O', 'O', 'C', 'C', 'C', 'H', 'H', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'C', 'O', 'C', 'O', 'C', 'H', 'H', 'O', 'H', 'H', 'C', 'O', 'O', 'O', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H']
Using float64 for MACECalculator, which is slower but more accurate. Recommended for geometry optimization.
/usr/local/lib/python3.10/dist-packages/mace/calculators/mace.py:135: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See <https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models> for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.
torch.load(f=model_path, map_location=device)
Default dtype float64 does not match model dtype float32, converting models to float64.
タンパク質の原子構造や各原子の座標などの情報が表示されています。
MDの実行
それでは今までに定義した関数を使って実際に分子動力学計算を実行してみましょう！
simpleMD(molecule, T_init=1200, calc=mace_calc, fname='mace_md.xyz', s=10, T=2000)
先ほど定義したsimpleMD関数を実行しています。引数は次の通りです。
molecule:- ASE の
Atomsオブジェクトで、分子の初期構造。 - ここでは、
read(pdb_path)で読み込まれたsolvent_molecs.xyzファイルの内容。
- ASE の
T_init=1200:- 初期温度（単位: ケルビン）。
- ここでは 1200 Kを設定。
calc=mace_calc:mace_calcは、mace_off関数で作成された計算設定。- 機械学習ベースの力場モデルを使用してエネルギーや力を計算。
fname='mace_md.xyz':- シミュレーション結果の出力ファイル名。
mace_md.xyzファイルとして保存。
s=10:- シミュレーション結果の出力間隔。
- 10ステップごとに出力データをファイルに保存。
T=2000:- 総シミュレーションステップ数。
- 2000ステップ分のMDを実行。
これを実行すると最終的に次のようなグラフが得られます。
このグラフを見ると、分子動力学計算実行中のエネルギーのブレは0.1eVほどであり、温度のブレも500Kほどです。単純な溶媒システム自体がMDシミュレーションで正常に動作していることを評価できました。お疲れ様でした。
最後に
MACEを行うことであなたのPCから簡単に分子動力学計算を実行でき、分子の動的特性や熱力学的挙動を理解することで化学や物理の問題に対する新しい知見を得ることができます。ぜひトライしてみてください！
参考文献
MACE tutorial https://github.com/ilyes319/mace-tutorials/blob/main/mace-users/MACE_users.ipynb
