Google Colab上でAtomic Simulation Environment(ASE)によるCuの分子動力学（MD）シミュレーション

このチュートリアルでは、ASE（Atomic Simulation Environment） を用いて、結晶構造の可視化・操作から、MDシミュレーションの実行・可視化・分析までを段階的に紹介しています。またGoogle Colab上でも実行可視化できる様にしています。

Atomic Simulation Environment (ASE)とは

第一原理計算/量子化学計算を行う超便利なPythonライブラリ
https://wiki.fysik.dtu.dk/ase/
計算の部分は既存のさまざまな計算パッケージを使用できる。今回はEffective medium theory (EMT)というポテンシャルの実装を用いる。

📦 使用ライブラリとインストール

最初に、必要なライブラリをインストールします。

!pip install ase pymatgen matplotlib py3Dmol ipywidgets

ase: シミュレーションの中心ライブラリ。
py3Dmol: 3D可視化ツール。
matplotlib: グラフ描画。
pymatgen: 結晶構造の解析に有用。

import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms, units
from ase.io import read, write
from ase.spacegroup import crystal
from ase.build import make_supercell, bulk
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.langevin import Langevin
from ase.constraints import FixAtoms
from ase.calculators.emt import EMT
from ase.io.trajectory import Trajectory
import py3Dmol
from google.colab import files
import os
import json

👀 py3Dmolを用いた原子構造の可視化

構造やTrajectoryの可視化に必要な関数を先に定義しておきます。
view_atoms_py3dmolはpy3Dmolを使用した構造描画用の関数です。

# py3Dmolを使った可視化関数
def view_atoms_py3dmol(atoms, width=500, height=400):
    """ASE原子オブジェクトをpy3Dmolで表示するためのヘルパー関数"""
    # 原子位置と元素シンボル、ユニットセル情報（あれば）を取得
    positions = atoms.get_positions()
    symbols = atoms.get_chemical_symbols()
    cell = atoms.get_cell()

    # XYZに変換
    xyz = f"{len(atoms)}\nAtoms\n"
    for i, (pos, sym) in enumerate(zip(positions, symbols)):
        xyz += f"{sym} {pos[0]:.6f} {pos[1]:.6f} {pos[2]:.6f}\n"

    # py3Dmolビューアを作成
    view = py3Dmol.view(width=width, height=height)
    view.addModel(xyz, "xyz")

    # 原子の表示スタイルを設定
    view.setStyle({'sphere':{'colorscheme':'Jmol', 'scale':0.3}})
    view.zoomTo()
    return view

# MDトラジェクトリを可視化するための関数
def view_trajectory_py3dmol(traj, start=0, end=None, step=1, width=500, height=400):
    """ASEトラジェクトリをpy3Dmolでアニメーション表示する関数"""
    if end is None:
        end = len(traj) #計算途中でも、正常に読み込めるように。

    # 最初のフレームを表示
    view = view_atoms_py3dmol(traj[start], width=width, height=height)

    # フレームをモデルとして追加
    for i in range(start+step, end, step):
        atoms = traj[i]
        # 原子位置と元素シンボルを取得
        positions = atoms.get_positions()
        symbols = atoms.get_chemical_symbols()

        # XYZフォーマットに変換して、append.
        xyz = f"{len(atoms)}\nFrame {i}\n"
        for j, (pos, sym) in enumerate(zip(positions, symbols)):
            xyz += f"{sym} {pos[0]:.6f} {pos[1]:.6f} {pos[2]:.6f}\n"

        view.addModel(xyz, "xyz")

    # すべてのモデルに同じスタイルを適用
    view.setStyle({'sphere':{'colorscheme':'Jmol', 'scale':0.3}})

    # アニメーション設定
    view.animate({'loop': 'forward', 'interval': 500})
    view.zoomTo()
    return view

なぜ .xyz を使うのか？
py3Dmol は様々なフォーマット（.pdb, .mol, .xyzなど）を扱えるものの、ASE の Atoms オブジェクトを直接的に扱えないので、一度文字列としてXYZフォーマットにする必要があります。

🧱 ステップ1: 結晶構造の取得と表示

銅（Cu）のFCC構造を定義。
.cif ファイルとして保存し、構造の基本情報（化学式、格子定数など）を表示。
3Dで構造も表示。

atoms = bulk('Cu', 'fcc', a=3.61) # aはLattice constant (格子定数)です。

# 構造をcifファイルとして保存
cif_filename = 'copper.cif'
write(cif_filename, atoms)
print(f"作成された構造: {cif_filename}")

# 原子構造の基本情報を表示
print("結晶構造の基本情報:")
print(f"化学式: {atoms.get_chemical_formula()}")
print(f"原子数: {len(atoms)}")
print(f"格子定数: {atoms.cell.cellpar()[:3]} Å")
print(f"格子角度: {atoms.cell.cellpar()[3:]} 度")

# py3Dmolを使用して3D可視化
view_atoms = view_atoms_py3dmol(atoms)
view_atoms.show()

作成された構造: copper.cif
結晶構造の基本情報:
化学式: Cu
原子数: 1
格子定数: [2.55265548 2.55265548 2.55265548] Å
格子角度: [60. 60. 60.] 度


画像1：py3dmolによる表示結果（Cu原子）

🧩 ステップ2: 結晶構造の操作

ここでは、Cu（1原子）だったものを周りに複製します。この様なものをsupercellと言います。

単位胞を3次元方向に拡張し、スーパーセルを作成し、全部で8原子の構造を作ります。ASEを使用すると、非常に簡単です。
原子数や格子のサイズを表示し、可視化します。

# スーパーセルの作成
nx = 2
ny = 2
nz = 2

# スーパーセルの作成（単位胞を指定方向に繰り返す）
supercell = atoms.repeat((nx, ny, nz))
print(f"スーパーセル情報:")
print(f"原子数: {len(supercell)}")
print(f"セルの大きさ: {supercell.cell.cellpar()[:3]} Å")

# スーパーセルの可視化
view_supercell = view_atoms_py3dmol(supercell)
view_supercell.show()

スーパーセル情報:
原子数: 8
セルの大きさ: [5.10531096 5.10531096 5.10531096] Å


画像1：py3dmolによる表示結果（Cu₈）

🔬 ステップ3: MDシミュレーションの実行

⚙️ 計算条件と初期設定

計算器には EMT（Effective Medium Theory）を使用。
温度: 500 K、圧力: 100 bar、時間ステップ: 5 fs。
インデックス0の原子を固定（FixAtoms）。
ランジュバン熱浴を使ったNVTアンサンブルでMD実行。NVTアンサンブルでは、原子数(N)、体積（V）、温度（T）が一定に保たれます。温度だけは、サンプリングするので、指定温度の近くで多少揺れます。
マクスウェル・ボルツマン分布から、初期速度として各原子の速度ベクトルをサンプリングします。

# MDのセットアップ
md_atoms = supercell.copy()
md_atoms.calc = EMT()

# 温度と圧力の設定
temperature = 500  # K
pressure = 100 * units.bar  # 100 bar

print(f"MDの設定:")
print(f"温度: {temperature} K")
print(f"圧力: {pressure / units.bar} bar")
print(f"原子数: {len(md_atoms)}")

# トラジェクトリファイルの設定
traj_filename = 'md_trajectory.traj'
traj = Trajectory(traj_filename, 'w', md_atoms)

# 初期速度を設定（マクスウェル・ボルツマン分布）
MaxwellBoltzmannDistribution(md_atoms, temperature_K=temperature)

# NVT分子動力学（ランジュバン動力学を使用）
time_step = 5.0  # fs
print(f"時間ステップ: {time_step} fs")

# MDインテグレータの設定
constraint = FixAtoms(indices=[0])  # インデックス0の原子を固定
md_atoms.set_constraint(constraint)
friction = 0.002  # 摩擦係数
dyn = Langevin(md_atoms, time_step * units.fs, temperature_K=temperature, friction=friction)

# データ保存用リスト
energies = []
temperatures = []
pressures = []

# データ保存関数
def store_data():
    energies.append(md_atoms.get_potential_energy())
    temperatures.append(md_atoms.get_temperature())
    pressures.append(pressure)  # ASEは圧力の直接取得が困難なため、設定値を保存

    # カスタムデータを traj に保存
    md_atoms.info['energy'] = energies[-1]
    md_atoms.info['temperature'] = temperatures[-1]
    md_atoms.info['pressure'] = pressures[-1]
    traj.write(md_atoms)

# トラジェクトリの保存（10ステップごと）
dyn.attach(store_data, interval=10)

# MDステップ数の設定
md_steps = 10000
print(f"MDステップ数: {md_steps}")
print("MDシミュレーションを実行中...")

# MDの実行
dyn.run(md_steps)
print("MDシミュレーション完了")

📊 データの記録と保存

10ステップごとに、温度や堆積の変化を.csv へ保存します。
.traj ファイルにStep数、全エネルギー、構造(原子座標)、速度ベクトルを保存します。中身は.xyzファイルの塊です。他にも電荷なども保存できます。

# 温度・エネルギー・圧力データを CSV に保存
data = pd.DataFrame({
    'Step': np.arange(0, len(energies) * 10, 10),  # 10ステップごと
    'Energy (eV)': energies,
    'Temperature (K)': temperatures,
    'Pressure (bar)': pressures
})

csv_filename = 'md_simulation_data.csv'
data.to_csv(csv_filename, index=False)
print(f"MDシミュレーションデータを保存しました: {csv_filename}")

# トラジェクトリをXYZ形式に変換
xyz_filename = 'md_trajectory.xyz'
write(xyz_filename, Trajectory(traj_filename), format='xyz')
print(f"トラジェクトリをXYZ形式に変換しました: {xyz_filename}")
print(f"トラジェクトリファイル {xyz_filename} をダウンロード")
files.download(xyz_filename) # Google Colabのprogram内からGoogle Driveにトラジェクトリファイルを保存

📈 ステップ4: MD結果の可視化

matplotlib を用いて以下の3つをプロットします。：
- エネルギー
- 温度
- 圧力


# MD ステップの範囲
md_steps_range = range(0, len(energies) * 10, 10)

# プロットの作成
plt.figure(figsize=(10, 8))

# 1. ポテンシャルエネルギー
plt.subplot(3, 1, 1)
plt.plot(md_steps_range, energies, label='Potential Energy (eV)', color='blue')
plt.xlabel('MD Step')
plt.ylabel('Energy (eV)')
plt.title('Potential Energy during MD Simulation')
plt.grid(True)
plt.legend()

# 2. 温度
plt.subplot(3, 1, 2)
plt.plot(md_steps_range, temperatures, label='Temperature (K)', color='red')
plt.xlabel('MD Step')
plt.ylabel('Temperature (K)')
plt.title('Temperature Fluctuation during MD Simulation')
plt.grid(True)
plt.legend()

# 3. 圧力
plt.subplot(3, 1, 3)
plt.plot(md_steps_range, pressures, label='Pressure (bar)', color='green')
plt.xlabel('MD Step')
plt.ylabel('Pressure (bar)')
plt.title('Pressure Variation during MD Simulation')
plt.grid(True)
plt.legend()

# グラフの調整と表示
plt.tight_layout()
plt.show()


画像3：エネルギー、温度、圧力の経時変化。エネルギー、温度は上昇していっていますね。温度は500Kあたりにありますが、まだ平衡化してるとは言えなさそうですね、10000 stepでは足りない様です。


最後にCu₈のシミュレーションの様子をgifにしたので見ていただければと思います。index=0の原子は固定してあります（真空中をふわふわしないように）。
このgifはOVITO(https://www.ovito.org/ )という可視化ツールを使用しています。OVITOの使い方解説は、またの機会に。

📝 おわりに

このスクリプトは、結晶構造の取得 → スーパーセルの構築 → MD実行 → 可視化 という一連の流れを体験できる、ASEのチュートリアルです。特にGoogle Colab上での実行を想定しているため、誰でも手軽に3D表示と解析を行うことができます。

化学徒としては、プログラム書くよりも、MDシミュレーションで化学構造やその変化を見る方が楽しいと思う。
でも、ある程度のプログラムは書けないと始まらないですけどね。。。