ASEで分子を作成

計算化学

PythonのASE (Atomic Simulation Environment) で分子を作成する際の基本について簡単にまとめたいと思います。

スポンサーリンク

<class ‘ase.atoms.Atoms’>

分子はこのクラスのオブジェクトです。試しに酢酸分子を作ってみましょう。

from ase.build import molecule

atoms = molecule("CH3COOH")
pos = atoms.get_positions()
elm = atoms.get_chemical_symbols()
num = atoms.get_atomic_numbers()
mas = atoms.get_masses()

print(atoms.symbols)
for i in range(atoms.get_global_number_of_atoms()):
    print("%3i %3s %3i %8.3f"%(i, elm[i], num[i], mas[i]), pos[i])

print(mas.sum())

この「atoms」が<class ‘ase.atoms.Atoms’>のオブジェクト、すなわち分子です。
このクラスで使えるメソッドはたくさんありますが、ここで使っているものは

  • .get_positions() 原子数 × XYZ座標の二次元配列
  • .get_chemical_symbols() 元素記号
  • .get_atomic_numbers() 原子番号
  • .get_masses() 原子量

出てくるものはNumpyArrayなので、例えば最後にprintしているように原子量の.get_masses()を.sum()すると分子量になります。

実行結果はこうなります。

CO2HCH3
  0   C   6   12.011 [0.      0.15456 0.     ]
  1   O   8   15.999 [0.166384 1.360084 0.      ]
  2   O   8   15.999 [-1.236449 -0.415036  0.      ]
  3   H   1    1.008 [-1.867646  0.333582  0.      ]
  4   C   6   12.011 [ 1.073776 -0.892748  0.      ]
  5   H   1    1.008 [ 2.048189 -0.408135  0.      ]
  6   H   1    1.008 [ 0.968661 -1.528353  0.881747]
  7   H   1    1.008 [ 0.968661 -1.528353 -0.881747]
60.05199999999999

なお、from ase.build import moleculemolecule("CH3COOH") のように名前で指定できるものの一覧はこちらです。

>>> from ase.collections import g2
>>> print(g2.names)
['PH3', 'P2', 'CH3CHO', 'H2COH', 'CS', 'OCHCHO', 'C3H9C', 'CH3COF',
 'CH3CH2OCH3', 'HCOOH', 'HCCl3', 'HOCl', 'H2', 'SH2', 'C2H2', 'C4H4NH',
 'CH3SCH3', 'SiH2_s3B1d', 'CH3SH', 'CH3CO', 'CO', 'ClF3', 'SiH4',
 'C2H6CHOH', 'CH2NHCH2', 'isobutene', 'HCO', 'bicyclobutane', 'LiF',
 'Si', 'C2H6', 'CN', 'ClNO', 'S', 'SiF4', 'H3CNH2', 'methylenecyclopropane',
 'CH3CH2OH', 'F', 'NaCl', 'CH3Cl', 'CH3SiH3', 'AlF3', 'C2H3', 'ClF',
 'PF3', 'PH2', 'CH3CN', 'cyclobutene', 'CH3ONO', 'SiH3', 'C3H6_D3h',
 'CO2', 'NO', 'trans-butane', 'H2CCHCl', 'LiH', 'NH2', 'CH',
 'CH2OCH2', 'C6H6', 'CH3CONH2', 'cyclobutane', 'H2CCHCN', 'butadiene',
 'C', 'H2CO', 'CH3COOH', 'HCF3', 'CH3S', 'CS2', 'SiH2_s1A1d', 'C4H4S',
 'N2H4', 'OH', 'CH3OCH3', 'C5H5N', 'H2O', 'HCl', 'CH2_s1A1d', 'CH3CH2SH',
 'CH3NO2', 'Cl', 'Be', 'BCl3', 'C4H4O', 'Al', 'CH3O', 'CH3OH', 'C3H7Cl',
 'isobutane', 'Na', 'CCl4', 'CH3CH2O', 'H2CCHF', 'C3H7', 'CH3', 'O3',
 'P', 'C2H4', 'NCCN', 'S2', 'AlCl3', 'SiCl4', 'SiO', 'C3H4_D2d',
 'H', 'COF2', '2-butyne', 'C2H5', 'BF3', 'N2O', 'F2O', 'SO2', 'H2CCl2',
 'CF3CN', 'HCN', 'C2H6NH', 'OCS', 'B', 'ClO', 'C3H8', 'HF',
 'O2', 'SO', 'NH', 'C2F4', 'NF3', 'CH2_s3B1d', 'CH3CH2Cl',
 'CH3COCl', 'NH3', 'C3H9N', 'CF4', 'C3H6_Cs', 'Si2H6', 'HCOOCH3',
 'O', 'CCH', 'N', 'Si2', 'C2H6SO', 'C5H8', 'H2CF2', 'Li2', 'CH2SCH2',
 'C2Cl4', 'C3H4_C3v', 'CH3COCH3', 'F2', 'CH4', 'SH', 'H2CCO', 
'CH3CH2NH2', 'Li', 'N2', 'Cl2', 'H2O2', 'Na2', 'BeH', 'C3H4_C2v', 'NO2']
スポンサーリンク

化合物名やSMILESやCIDで分子を作成

PubChemを使って分子を作ることもできます。

from ase.data.pubchem import pubchem_atoms_search

benzene1 = pubchem_atoms_search(name="benzene")
benzene2 = pubchem_atoms_search(cid=241)
benzene3 = pubchem_atoms_search(smiles="c1ccccc1")
print(benzene1.symbols)
print(benzene1 == benzene2, benzene1 == benzene3)

結果はこうなります。

C6H6
True True

化合物名、CID、SMILES、いずれも同じオブジェクトが作られました。

また、複数の配座がある場合は、

from ase.data.pubchem import pubchem_atoms_search, pubchem_atoms_conformer_search

octane_conformers = pubchem_atoms_conformer_search(name="octane")

このoctane_conformersはPythonのリストなので、octane_conformers[0]のようにして1個目の分子のオブジェクトを選択できます。

スポンサーリンク

座標を入力して作る場合

二酸化炭素の分子を作ってみます。

from ase import Atoms

CO2 = Atoms("OCO", positions=[(-1.16, 0, 0), (0, 0, 0), (1.16, 0, 0)])
print(CO2.get_chemical_symbols())
print(CO2.get_masses().sum())
print(CO2.get_all_distances())

結果はこうなります。

['O', 'C', 'O']
44.009
[[0.   1.16 2.32]
 [1.16 0.   1.16]
 [2.32 1.16 0.  ]]

最後の.get_all_distances()は各原子同士の距離を出力します。
2番目の[1.16 0. 1.16]を見ると、Cの原子が2つのOに対して1.16Åの距離にあることが分かります。

スポンサーリンク

可視化

先ほどのCO2分子を可視化します。

from ase import Atoms
from ase.visualize import view

CO2 = Atoms("OCO", positions=[(-1.16, 0, 0), (0, 0, 0), (1.16, 0, 0)])
view(CO2)

実行するとViewerが出てきます。

View設定をするとこんな感じ。

(Bond情報ってどこにあるんだろう…)

スポンサーリンク

ファイルのインポート

GAMESS_USの計算結果をインポートしてみます。

from ase.io import read

atoms = read("benzene.out")
print(atoms.get_atom_positions())

結果はこうなります。

[[ 0.64489995 -1.23569991  0.        ]
 [ 1.3927999  -0.0594      0.        ]
 [ 0.74779995  1.17629991  0.        ]
 [-0.64489995  1.23579991  0.        ]
 [-1.3926999   0.0595      0.        ]
 [-0.74789995 -1.17639991  0.        ]
 [ 1.14799992 -2.20039984  0.        ]
 [ 2.47939982 -0.10579999  0.        ]
 [ 1.3313999   2.09439985  0.        ]
 [-1.14829992  2.20029984  0.        ]
 [-2.47959982  0.10589999  0.        ]
 [-1.3314999  -2.09449985  0.        ]]

このoutファイルはGAMESSの構造最適化の結果なのですが、実は最適化された座標ではなく初期座標が読み込まれています。
構造最適化の結果の座標を読み込むことはできないようです。
GAMESSで構造最適化をしたい場合はASE上で構造最適化を実行する必要があります。
詳しくはこちらの記事で解説しています。

コメント

タイトルとURLをコピーしました