Colabで分子モデリングからMLPで計算【MOF編】

0. セットアップ（最初に一度だけ）

これから UiO-66 の構造に CO₂ を入れていきます。

でもその前に、“道具箱”を整えるところから始めましょう。

たとえば実験なら、フラスコや試薬を並べてから反応を始めますよね。

同じように、ここでは分子を作ったり並べたりするための Python ツールをインストールしておきます。

用意するライブラリたち

• ASE：構造を読み書きしたり、結晶のセルを扱う“万能助手”。

  CIF の読み書きや、後で UMA に渡すのもこの子が担当です。

• RDKit：SMILES から分子の3D構造をつくる、化学屋の味方。

  「O=C=O」と打つだけで CO₂ の立体構造を出してくれます。

• pymatgen：CIF の秩序や占有率をチェックする監査官。

  実験CIFにありがちな「0.5サイト混在」も見逃しません。

セットアップコード

まずはこれを一度だけ実行すれば準備完了です。

!pip -q install ase rdkit pymatgen # Colabで失敗する場合は: pip install rdkit-pypi

import numpy as np, os, pathlib
from ase.io import read, write
from ase import Atoms
from rdkit import Chem
from rdkit.Chem import AllChem
# RDKit→ASE（単一コンフォーマ）変換
def smiles_to_rdkit_3d(smiles, num_conf=10, max_iters=200): mol = Chem.AddHs(Chem.MolFromSmiles(smiles)) params = AllChem.ETKDGv3(); params.randomSeed = 42 AllChem.EmbedMultipleConfs(mol, numConfs=num_conf, params=params) res = AllChem.MMFFOptimizeMoleculeConfs(mol, maxIters=max_iters) best_id = int(np.argmin([r[1] for r in res])) conf = mol.GetConformer(best_id) m = Chem.Mol(mol); m.RemoveAllConformers(); m.AddConformer(conf, assignId=True) return m
def rdkit_to_ase(mol): conf = mol.GetConformer() syms, pos = [], [] for a in mol.GetAtoms(): p = conf.GetAtomPosition(a.GetIdx()) pos.append([p.x, p.y, p.z]); syms.append(a.GetSymbol()) return Atoms(symbols=syms, positions=np.array(pos, float))
def to_large_p1_cell(atoms, box=30.0, pbc=True): a = atoms.copy(); a.set_cell(); a.set_pbc((pbc,pbc,pbc)); a.center() return a

この関数たちがやっていること

• smiles_to_rdkit_3d

  → SMILES 文字列をもとに 複数の立体候補を生成し、最も安定な構造を選ぶ。

  乱数シードを固定しているので、何度やっても同じ結果になります。

• rdkit_to_ase

  → RDKit 分子を ASE の Atoms に変換。これで「実験室」から「計算室」に渡せるようになります。

• to_large_p1_cell

  → 分子を P1 の立方セルにポンと置く。周期境界をつけておくと、後で UMA や DFT コードにそのまま投げられます。

これで準備完了です。

次のステップでは、UiO-66（MOF骨格）を読み込み、CO₂ の3D構造を作って並べてみる

ところへ進みましょう。いよいよ“分子を置く”作業が始まります。

1. 初期構造を生成 / 読み込み（UiO-66 を持ってくる）

1-1. CIF を Colab に置く

方法A：手でアップロード（確実）

1. 画像どおりに COD 検索（https://www.crystallography.net/cod/search.php）で UiO-66 を入れる
2. 一覧から 4132636.cif をクリックして保存
3. Colab 左のフォルダアイコン → Upload → 4132636.cif を選ぶ

方法B：コマンドで取得（速い）

多くの COD エントリは https://www.crystallography.net/cod/<ID>.cif で直接落とせます。

もし 404 になる場合は方法Aで保存してください。

!wget -O 4132636.cif https://www.crystallography.net/cod/4132636.cif || echo "直接DLできない場合は手動へ"

小ワザ：作業しやすいように、以降は mof.cif という共通名で扱います。

import shutil, os
if os.path.exists("4132636.cif"): shutil.copy("4132636.cif", "mof.cif")
print("exists:", os.path.exists("mof.cif"))

1-2. まずは“中身をチラ見”して安心する

CIF はテキストです。先頭だけ眺めると、由来やセルがすぐわかる。

!head -n 45 mof.cif

ここで 化学式・格子定数・空間群 がそれっぽいか軽く確認。UiO-66 はしばしば Fm-3m（No.225） が出ます。

1-3. ASE/pymatgenで読み込み＆基本情報

最初の読み込みで警告が出ても、単点差分が目的なら気にしすぎなくてOK。

from ase.io import read
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
atoms = read("mof.cif") # 周期系として読む
structure = Structure.from_file("mof.cif")
print("Natoms:", len(atoms))
print("Cell lengths (Å):", atoms.cell.lengths())
print("Cell angles (deg):", atoms.cell.angles())
sga = SpacegroupAnalyzer(structure, symprec=1e-2)
print("Space group:", sga.get_space_group_symbol(), sga.get_space_group_number())

もし 占有率 < 1 の“混在サイト”が含まれていたら（実験CIFあるある）、後で is_ordered を見ます。単一点評価なら P1 に落として使えば実務的には十分です。

1-4. その場でサッと可視化

可視化ツール を立ち上げなくても、雰囲気はこれで掴めます。

from ase.io import read, write
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
# 人が見やすい2Dスナップを作成
mof_plus_guest = read("mof.cif")
fig, ax = plt.subplots(figsize=(5, 4.5), dpi=250)
# 回転角度: 'degXx,degYy,degZz' 形式（例：'0x,45y,0z'）
# 原子半径: radii=0.9 前後で調整（0.3〜1.2くらいが目安）
plot_atoms( mof_plus_guest, ax, rotation='0x,0y,0z', radii=0.9, show_unit_cell=2 # セル枠も薄く表示（=2で枠のみ）
)
ax.set_axis_off()
plt.tight_layout()

もっと拡大して見たい／原子を消したい時は、のちほど Avogadro でも確認します。

1-5. “作業用プレビュー”を書き出す（後工程の定番）

以降の章で編集しても原本を壊さないよう、コピーを作ってそちらを触ります。

from ase.io import write
framework_preview = "framework_preview.cif"
write(framework_preview, atoms)
print("wrote:", framework_preview)

この framework_preview.cif を使って、CO₂ を入れる・回す・位置を変えるといった作業に進みます。

（※原本 mof.cif はそのまま温存しておくと、いつでもやり直せて安心。）

1-6. 補足：COD以外で MOF CIF を探すなら？

💬 個人的おすすめは CoRE-MOF（“Computation-Ready Experimental MOF Database”）。

UiO-66 を含む主要MOFが、溶媒除去済みでそのまま計算に使える からです。

ここまでで、UiO-66 の CIF が Colab に入り、読む・見る・複製するの三拍子がそろいました。CO₂ の 3D を作って、いよいよ“置いてみる”に進みます。

2. 軽い微調整：Avogadro で“整える”

ここで、MOFに入れるゲスト分子を作ります。今回は下記のようにCO2を作成します。

from ase.build import molecule
from ase.io import write
# CO2分子を作成し、guest_co2.xyzとして保存
co2 = molecule("CO2")
write("guest_co2.xyz", co2)
print("CO2ファイルを作成しました → guest_co2.xyz")

他の分子を作成したい場合、下記のコードで作成できる分子を確認して下さい。

かなりの分子をASEを使って作成することが出来ます。

from ase.collections import g2
print(g2.names) #doctest:+NORMALIZE_WHITESPACE

次に作成したファイルをAvogadro or VESTAで確認していきます。

Avogadroは以下の公式サイトでダウンロードすることが出来ます。

Avogadroのダウンロードページを開くと以下のようなページに遷移するので、タブ内中央の”Download”をクリックしてください。自動的にダウンロードが開始されます。

Avogadroの詳細なダウンロード方法・使い方は下記の記事を参考にしてください。。

【Avogadro】AvogadroでGaussianの構造最適化入力ファイルを作成しよう！ – LabCode

VESTAは下記のダウンロードサイトからセットアップしてください。

Download – VESTA

ここは“髪を整える”くらいの小さな手直し。触りすぎないのがコツ。

分子（CO₂）

1. guest_co2.xyz を Avogadro で開く
2. 直線性と衝突を目視
3. 必要なら Extensions → Optimize Geometry → UFF or MMFF を軽く一回
4. そのまま上書き保存（ファイル名は維持）

骨格（UiO-66）

1. framework_preview.cif を Avogadro か VESTA で開く
2. 格子と配位環境を眺めるだけ（空間群は触らない）
3. 保存する場合は CIF形式のまま。ここでも“触りすぎない”

❗ 骨格の μ₃-OH などは構造の一部なので削らないでOK。溶媒だけを落とす場合は、元から溶媒なし構造を取り直す方が安全です。

P1運用（必要になったら）

# 対称性で“跳ぶ”挙動が気になるときは P1 に落とす
from ase.io import read, write
tmp = read("framework_preview.cif")
write("framework_P1.cif", tmp) # 実質 P1 相当で出力

VESTAで可視化したUiO-66

3. 置いてみる：ゲスト配置 → 合成CIFをつくる

いよいよ“分子を置く”時間。分率座標で吸着サイトを狙い撃ち。

# MOF骨格のプレビューCIF（必要ならファイル名を整えるだけ）
framework_preview = "framework_preview.cif"
# Load the guest molecule from the previously created xyz file
guest = read("guest_co2.xyz")
# まずはセル中心付近に1分子置く簡単例
combined = read(framework_preview).copy()
cell = combined.get_cell()
center_cart = (cell[0] + cell[1] + cell[2]) * 0.5
mol = guest.copy()
for ax in ("x","y","z"): # ランダム回転（必要なら固定角度に） if ax == "x": axis = (1, 0, 0) elif ax == "y": axis = (0, 1, 0) elif ax == "z": axis = (0, 0, 1) mol.rotate(np.random.uniform(0,360), v=axis, center="COM")
mol.translate(center_cart - mol.get_center_of_mass())
combined += mol
out_cif = "mof_plus_guest.cif"
write(out_cif, combined)
print("wrote:", out_cif)
# 分率座標指定で配置したいときの関数
def place_guest_by_frac(framework_atoms, guest_atoms, frac_xyz=(0.5,0.5,0.5)): f = framework_atoms.copy() g = guest_atoms.copy() # 分率→直交 frac = np.array(frac_xyz) cart = frac @ f.get_cell().array # (1x3) * (3x3) # 回転など必要ならここで g.translate(cart - g.get_center_of_mass()) return f + g
# 例：a/4, b/4, c/2 に置く
combined2 = place_guest_by_frac(read(framework_preview), guest, (0.25,0.25,0.25))
write("uio66_co2_safe.cif", combined2)

今回は分率座標指定で作成しています。CO2をa/4, b/4, c/4 に置いて、フレームワークと合体したファイルをuio66_co2_safe.cifというファイルで出力しています。

下記で分子が入っているか確認してみましょう。

4. 目視スナップ（2D図を作って保存）

下記を実行するとCO2がMOFのフレームワーク内に入っていることがわかります。

分かりにくい場合は、回転させて見やすい位置を探してください。

# 必要なimport（未済なら）
from ase.io import read, write
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
# ① P1に落としておく（分率ズレや対称性の自動補正を避ける）
# ※「分率座標で狙う」で書き出したファイル名に合わせて置き換えてOK
tmp = read("uio66_co2_safe.cif") # ← 直前に作った合成CIF
write("mof_plus_guest_frac_P1.cif", tmp) # ← P1運用ファイル
# ② 人が見やすい2Dスナップを作成
mof_plus_guest = read("mof_plus_guest_frac_P1.cif")
fig, ax = plt.subplots(figsize=(5, 4.5), dpi=250)
# 回転角度: 'degXx,degYy,degZz' 形式（例：'0x,45y,0z'）
# 原子半径: radii=0.9 前後で調整（0.3〜1.2くらいが目安）
plot_atoms( mof_plus_guest, ax, rotation='0x,45y,0z', radii=0.9, show_unit_cell=2 # セル枠も薄く表示（=2で枠のみ）
)
ax.set_axis_off()
plt.tight_layout()
plt.savefig("uio66_co2_view.png", bbox_inches="tight", pad_inches=0.02)
print("saved: uio66_co2_view.png")

コツ

• 回転は 'ax,by,cz' の角度指定。吸着サイトが見える方向を探す感覚で何枚か出すと◎

  例：'30x,30y,0z', '0x,60y,0z', '45x,45y,30z'

• radii は印象が変わりやすいので、0.7–1.0あたりを中心に微調整。

• show_unit_cell=2 でセル枠が入り、MOFのポア感が伝わりやすくなります。

• 複数候補を自動出力するなら、回転角だけforループで回して

  uio66_co2_view_01.png, _02.png, ... として保存しておくと便利 (付録で記載)。

5. 最終チェック：距離・占有率・セル

ここは“出荷前の検査”。衝突がないか／占有率が健全か／セルが変になってないかを、数十秒で確かめます。

from ase.io import read
import numpy as np
from pymatgen.core import Structure
final_path = "mof_plus_guest_frac_P1.cif" # 確認したい合成CIF
atoms = read(final_path)
# 5-1) 最近接距離（PBC, 最小鏡像）— 衝突の赤信号を消す
D = atoms.get_all_distances(mic=True)
np.fill_diagonal(D, 1e9)
dmin = float(D.min())
print(f"[Min distance] {dmin:.3f} Å ← ~0.61 Å 未満なら衝突疑い・配置を見直す")
# 5-2) セルの基本情報 — “変なセル”になっていないか
a, b, c = atoms.cell.lengths()
alpha, beta, gamma = atoms.cell.angles()
print(f"[Cell lengths] a={a:.3f}, b={b:.3f}, c={c:.3f} Å")
print(f"[Cell angles ] α={alpha:.2f}, β={beta:.2f}, γ={gamma:.2f}°")
# 5-3) 占有率チェック — 実験CIF由来の“半占有”に注意
s = Structure.from_file(final_path)
print("[is_ordered] True=すべて占有率1.0:", s.is_ordered)
if not s.is_ordered: print(" ↳ 注意：無秩序サイトあり（占有率≠1.0）。単一点ならP1運用で実務的にはOK、厳密には展開を検討。") for i, site in enumerate(s.sites): if not site.is_ordered: print(f" - site #{i}: {site.species}") # 例: {'O': 0.5, 'N': 0.5} など

読み取り方（サッと）

• 最短距離が ~0.61 Å 未満 → どこかがめり込んでる合図。回転・初期高さを少し変えて再出力してください(しきい値 0.61 Å はプレーンのUiO-66の最短距離）。
• is_ordered = False → 実験CIFあるあるの“半占有”。P1固定で単点を回すだけなら許容だが、厳密な解析では明示展開（どちらかに寄せる）を。
• セルは UiO-66 らしい立方に近い数値が出ていればOK（細かな差は測定由来のことも多い）。

付録：UMA で単点（E_ads）— “数を出して感触を掴む”

定義（eV）

Eads=E[MOF+CO2]−E[MOF]−E[CO2]

負 → 吸着“有利”（発熱）／ 正 → 吸着“不利”（吸熱）

準備

下記のコードを実行してください。

機械学習ポテンシャルUMAの使用許可の取得やトークン作成は詳細は下記の記事を参考にしてください。

MOF吸着エネルギー計算の新手法：DFT精度を瞬時に再現【機械学習ポテンシャル】 – LabCode

!pip install fairchem-core

!huggingface-cli login

A. 最小例（1配置でサクッと試す）

from fairchem.core import FAIRChemCalculator, pretrained_mlip
from ase.io import read
# UMA の読み込み（ODAC: CO2/H2O 吸着向けタスク）
unit = pretrained_mlip.get_predict_unit("uma-s-1p1")
calc = FAIRChemCalculator(unit, task_name="odac") # device は自動。CPU固定なら device="cpu"
# 構造
mof_co2 = read("mof_plus_guest_frac_P1.cif") # UiO-66 + CO₂（P1）
mof_empty = read("mof.cif") # UiO-66（空）
co2 = read("guest_co2.xyz") # CO₂（分子）
for a in (mof_co2, mof_empty, co2): a.calc = calc
E_mof_co2 = mof_co2.get_potential_energy()
E_mof_empty = mof_empty.get_potential_energy()
E_co2 = co2.get_potential_energy()
Eads = E_mof_co2 - E_mof_empty - E_co2
print(f"Eads = {Eads:.3f} eV (negative=adsorption favorable)")

B. 複数候補から“まっとうな”ものを選ぶ（おすすめ）

配置次第で符号すら変わるのが実務。まずは複数配置を作り、最短距離で弾く→Eads最小を採用するとハズレにくい。

import os, numpy as np, itertools, csv, glob
from pymatgen.core import Structure, Molecule, Site
from pymatgen.core.operations import SymmOp
# ---------- 既存の関数（そのまま利用） ----------
def place_molecule_with_clearance(framework: Structure, mol: Molecule, frac_xyz=(0.5,0.5,0.5), min_clear=2.0, tries=300): lat = framework.lattice target_cart = lat.get_cartesian_coords(frac_xyz) for _ in range(tries): m = mol.copy() # ランダム回転 axis = np.random.rand(3); axis /= np.linalg.norm(axis) angle = np.random.rand()*360 m.apply_operation(SymmOp.from_axis_angle_and_translation(axis, angle, [0,0,0])) # COMを目標位置へ平行移動 shift = target_cart - m.center_of_mass m.translate_sites(range(len(m)), shift) # 一時的に合体して最近接をPBCでチェック tmp = framework.copy() start_idx = len(tmp) for sp, xyz in zip(m.species, m.cart_coords): tmp.append(sp, tmp.lattice.get_fractional_coords(xyz)) # PBC-aware ok = True for i in range(start_idx, len(tmp)): # ゲスト原子の近傍にフレーム原子が min_clear 未満で居ないか nns = tmp.get_neighbors(tmp[i], r=min_clear) if any(nn.index < start_idx for nn in nns): ok = False; break if ok: return tmp # 衝突なし → 採用 raise RuntimeError("配置に失敗（min_clearを下げる/位置を変える/triesを増やす）")
# ------------------------------------------------
# 入力（溶媒なし UiO-66 を推奨：CoRE-MOF / Materials Project / QMOF 等から）
framework_path = "/content/4132636.cif" # ここはあなたのUiO-66に合わせて
guest_path = "/content/guest_co2.xyz" # 3章で用意したCO2
mof = Structure.from_file(framework_path)
co2 = Molecule.from_file(guest_path)
# 置いてみる候補の“分率座標”たち（必要に応じて増やす/変更OK）
frac_points = [ (0.50, 0.50, 0.50), (0.25, 0.25, 0.25), (0.75, 0.25, 0.50), (0.25, 0.10, 0.25), (0.40, 0.40, 0.60), (0.60, 0.40, 0.40),
]
min_clear = 0.61 # Å：衝突回避の最小距離
tries = 300 # 試行回数（当たりが少ない座標は増やす）
os.makedirs("/content/candidates", exist_ok=True)
saved = []
for i, frac in enumerate(frac_points, start=1): try: placed = place_molecule_with_clearance(mof, co2, frac_xyz=frac, min_clear=min_clear, tries=tries) out = f"/content/candidates/uio66_co2_safe_{i:02d}.cif" placed.to(filename=out) print("saved:", out, "at frac", frac) saved.append(out) except RuntimeError as e: print(f"[skip] frac={frac}: {e}")
print("----")
print("safe candidates:", len(saved), "files")

作成したファイルを計算します。

import glob, csv
import numpy as np
from ase.io import read
from fairchem.core import FAIRChemCalculator, pretrained_mlip
# UMAの読み込み
unit = pretrained_mlip.get_predict_unit("uma-s-1p1")
calc = FAIRChemCalculator(unit, task_name="odac") # task_name=odac: gas adsorption向け
# 参照系（空MOFとCO2）
mof_empty = read("/content/4132636.cif")
co2 = read("/content/guest_co2.xyz")
for a in (mof_empty, co2): a.calc = calc
E_mof_empty = mof_empty.get_potential_energy()
E_co2 = co2.get_potential_energy()
print(f"E_mof_empty = {E_mof_empty:.3f} eV")
print(f"E_co2 = {E_co2:.3f} eV")
# 候補フォルダ内のCIFを探索
paths = sorted(glob.glob("/content/candidates/uio66_co2_safe_*.cif"))
print("Found", len(paths), "candidates.")
rows = []
for path in paths: a = read(path) a.calc = calc E_mof_co2 = a.get_potential_energy() Eads = E_mof_co2 - E_mof_empty - E_co2 rows.append((path, Eads)) print(f"{path.split('/')[-1]:30s} Eads = {Eads:7.3f} eV")
# 吸着エネルギーで並べ替え
rows.sort(key=lambda x: x[1])
# 保存
with open("/content/candidates/Eads_scan.csv", "w", newline="") as f: w = csv.writer(f) w.writerow(["path","Eads_eV"]) w.writerows(rows)
print("→ saved: /content/candidates/Eads_scan.csv")
# 最も安定な候補
best_path, best_E = rows[0]
print(f"\n✅ 最も安定: {best_path} (Eads = {best_E:.3f} eV)")

出力

• uio66_co2_safe_01.cif Eads = -1.399 eV
• uio66_co2_safe_02.cif Eads = -1.551 eV
• uio66_co2_safe_03.cif Eads = 92.727 eV
• uio66_co2_safe_04.cif Eads = 1.492 eV
• uio66_co2_safe_05.cif Eads = -1.307 eV
• uio66_co2_safe_06.cif Eads = -1.310 eV`

どう読む？

• 初期（“デフォルト”配置） → Eads ≈ -1.399 eV（非常に強い吸着）

• 位置調整後 ex. uio66_co2_safe_03.cif → Eads ≈ +92.727 eV（明確に不利）

  • 恐らくめり込んでいる。

• つまり 配置（姿勢・高さ）の違いだけで、符号が反転することがあります。

  → “1配置で結論を出さない” が鉄則。複数候補＋（必要なら）ごく軽い緩和で最小Eadsを採用するのが現実解。

目安：UiO-66×CO₂なら、文献や既知の“妥当帯域”は –0.3 〜 –0.6 eV 付近に落ちることが多いです。

−1 eV 級の強い値や +1 eV の正値は、初期配置のバイアスや局所衝突/反発の影響を疑って、構造最適化行い、配置や姿勢を見直してください。

最後に

いかがでしたか？本記事で紹介した分子構造の作り方を一度試していただき、機械学習ポテンシャルを使って計算を回していただくと、計算化学がより身近に感じられると思います。Cifファイルは論文などによく記載されているので、ゲスト分子を入れずに計算するだけでも面白いと思いますので是非試してみてください。

参考文献

1. K. Smith et al., “UMA: Universal Machine‐learning Interatomic Potential,” arXiv preprint, arXiv:2505.08762 (2025). https://arxiv.org/abs/2505.08762
2. FAIRchem Team, “FAIRChem: FAIR Chemistry Toolkit,” GitHub Repository (2025). https://github.com/facebookresearch/fairchem
3. Meta FAIR Chemistry, “facebook/UMA,” Hugging Face Model Card (2025). https://huggingface.co/facebook/UMA
4. B. M. Wood et al., “UMA: A Family of Universal Models for Atoms,” arXiv preprint, arXiv:2506.23971 (2025). https://arxiv.org/abs/2506.23971
5. S. Kitagawa, R. Kitaura, and S. Noro, “Functional Porous Coordination Polymers,” Angew. Chem. Int. Ed. 43, 2334–2375 (2004). https://doi.org/10.1002/anie.200300610
6. S. Horike, S. Shimomura, and S. Kitagawa, “Soft porous crystals,” Nature Chem. 1, 695–704 (2009). https://doi.org/10.1038/nchem.444
7. H. Furukawa, K. E. Cordova, M. O’Keeffe, and O. M. Yaghi, “The chemistry and applications of metal–organic frameworks,” Science 341, 1230444 (2013). https://doi.org/10.1126/science.1230444
8. J. H. Lee et al., “Design strategies for metal–organic frameworks in direct air capture of CO₂,” Coord. Chem. Rev. 463, 214536 (2022). https://doi.org/10.1016/j.ccr.2022.214536
9. J. H. Cavka et al., “A new zirconium inorganic building brick forming metal–organic frameworks with exceptional stability,” J. Am. Chem. Soc. 130, 13850–13851 (2008). https://doi.org/10.1021/ja8057953 PubMed
10. Crystallography Open Database (COD), Entry 4132636: UiO-66 (accessed 2025-11-03). https://www.crystallography.net/cod/4132636.cif crystallography.net