Opt¶
この章では系の局所安定座標とそのエネルギーを計算する、 構造最適化 (構造緩和などとも呼びます) について学びます。
最適化アルゴリズムには大域的最適化 (Global optimization)と、局所最適化(Local optimization)がありますが、構造最適化では局所最適化アルゴリズムを用います。
ここでは横軸は系の様々な配座(Conformation) を表し、縦軸はエネルギーを表します。
分子や結晶などの構造最適化を考える場合、N個の原子からなる系の座標値は 3*N 次元の自由度があり、例えば100原子からなる系では300次元もの高次元空間での最適化問題となります。 大域最適解を見つけることは一般に難しい問題です。
構造最適化は、系の初期座標を定めて、そこから局所安定点を見つけることにより、その安定構造を見つけます。
ある物質の組成(元素の組み合わせとその割合、例えばNaCl など)がわかっている場合でも、その構造がどの様になっているかは自明ではありませんが、 物質の多くの性質が安定構造付近の状態から決まるため、安定構造を求めることは重要です。
実験的に構造を決定する場合は以下のような分析機器・計測装置を使う事が考えられます。 しかし、その場合でも
スペクトル情報から物質の部分的な情報が得られるのみで、3次元構造の情報が直接観測できるわけではない (1~6)
単結晶や、薄膜といった構造を作成する必要があり、現実ではそういった試料の準備が難しいことがある (3, 8)
原子レベルのスケールで見るには解像度が足りない (7)
表面の形状は観察できるものの、その原子・分子の種類を区別することができない (9)
など、課題があり、今扱っている物質の各原子がどのように配置されているか、3次元構造の特定は容易ではありません。
分析機器の例：
名前 |
特徴 |
|---|---|
|
分子の振動による赤外線吸収を測定することで、分子構造の情報を得る |
|
ラマン散乱光を用いて物質の評価を行う |
|
結晶構造を持つ物質にX線を照射することで得られる回折線を解析することで物質の同定や定量を行う |
|
軟X線を物質に照射し、物質のイオン化に伴い放出される光電子e-を補足しエネルギー分析を行う |
|
物質を構成している個々の化合物の質量を正確に測定 |
|
原子核の共鳴現象を利用して、有機化合物の分析を行う |
|
波長の短い電子線を利用して0.4~4nmの解像度で表面を観察 |
|
電子線を当て透過した電子線の強弱により微小な物質を観察 |
|
探針と試料に作用する原子間力を検出する顕微鏡 |
分析機器の参考文献
そのため、シミュレーション上でどの構造になっていそうかを推定し、実験で検証する（またはその逆）を行うことが重要となります。
Opt molecules¶
まずは簡単な例として、周期境界条件のない単体の有機分子の構造最適化を行ってみます。
1. 水素分子¶
まずは、前章でも行った水素分子の最安定構造の結合距離を、構造最適化を用いて求めてみます。
以下のような流れとなります。
Step 1: ポテンシャルエネルギーを計算するCalculatorを用意
Step 2: 構造最適化を行いたい分子を用意
Step 3: 構造最適化アルゴリズムの手法を選択し、実行
[1]:
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
print(f"pfp_api_client: {pfp_api_client.__version__}")
estimator = Estimator(calc_mode=EstimatorCalcMode.MOLECULE, model_version="v3.0.0")
#estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v3.0.0")
calculator = ASECalculator(estimator)
pfp_api_client: 1.5.0
※PFP Estimator の calc_mode は以下のように用います。
一般の系：
EstimatorCalcMode.CRYSTAL_U0Uなしモード酸化物など：
EstimatorCalcMode.CRYSTALUありモード単体有機分子：
EstimatorCalcMode.MOLECULE分子モード
Step 2: 構造最適化を行いたい分子を用意
ここでは水素分子を作成します。座標値は人手で、1.0A 離れたところで初期座標を設定してみます。
[2]:
from ase import Atoms
atoms = Atoms("H2", [[0, 0, 0], [0, 0, 1.0]])
atoms.calc = calculator
Step 3: 構造最適化アルゴリズムの手法を選択し、実行
局所最適化アルゴリズムとして、ASEではNewton法に似た振る舞いで最適化を行う FIRE, MDMin や、2次関数近似をおこなって最適化を行う BFGS, LBFGS, BFGSLineSearch=QuasiNewton, LBFGSLineSearch法などが提供されています。これらのアルゴリズムの違いや性質については次節で説明をします。
今回は、LBFGS を用いて構造最適化を行ってみます。 基本的には atoms を局所最適化アルゴリズムに設定し、 run関数を呼ぶと構造最適化が行われます。 trajectory 引数を指定しておくと、構造最適化の過程を保存できます。
[3]:
import os
from ase.optimize import LBFGS, BFGS, FIRE
os.makedirs("output", exist_ok=True)
opt = LBFGS(atoms, trajectory="output/bfgs_h2.traj")
opt.run(fmax=0.01)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 02:19:43 -3.890512* 4.4993
LBFGS: 1 02:19:43 -4.402221* 3.2330
LBFGS: 2 02:19:43 -3.401438* 16.5629
LBFGS: 3 02:19:43 -4.551079* 2.2355
LBFGS: 4 02:19:43 -4.610468* 1.3879
LBFGS: 5 02:19:43 -4.641043* 0.3552
LBFGS: 6 02:19:43 -4.642630* 0.0576
LBFGS: 7 02:19:43 -4.642674* 0.0023
[3]:
True
今回の簡単な例では 7step、Matlantis では約2秒で構造最適化が終了しました。
構造最適化途中の系の可視化を行ってみます。
[4]:
from ase.visualize import view
from ase.io import Trajectory
from pfcc_extras.visualize.view import view_ngl
traj = Trajectory("output/bfgs_h2.traj")
view_ngl(traj, representations=["ball+stick"], w=400, h=300)
構造最適化により、原子距離がどのように推移したかをプロットしてみます。
[5]:
import numpy as np
# np.linalg.norm?
distance_list = []
E_pot_list = []
for atoms in traj:
p1, p2 = atoms.get_positions()
distance = np.linalg.norm(p1 - p2)
distance_list.append(distance)
E_pot_list.append(atoms.get_potential_energy())
print("Distance trajectory", distance_list)
Distance trajectory [1.0, 0.8714471953255789, 0.5432643813289825, 0.8178499966054394, 0.7851965174194478, 0.7317306649684048, 0.7426249533569255, 0.7411042704555475]
[6]:
import matplotlib.pyplot as plt
plt.plot(distance_list)
plt.xlabel("step")
plt.ylabel("Distance [A]")
plt.title("H2 opt trajectory")
plt.show()
最終的に、 0.74 A あたりで落ち着くことが確認できました。
前章ではすべての結合距離に対するエネルギーを計算し、最低エネルギーになる値を求めましたが、構造最適化を用いることで効率よく安定構造を見つけることができます。
2. シクロヘキサン¶
次の例として、シクロヘキサンの安定構造を考えてみます。
シクロヘキサンはC6H12で表される構造ですが、いくつかの立体配座(Conformation) を取ることが知られています。 代表的な配座異性体としては、以下の図のいす形 (A)、ねじれ舟形 (B)、舟形 (C)、半いす形 (D)などがあります。
一番安定な（エネルギーの低い）形はいす形で、25℃ではシクロヘキサン溶液中の全分子の99.99%がいす形になります。 ここでは、いす形 (A)および、舟形 (C)の初期構造を用意して、構造最適化を行い、その形やポテンシャルエネルギーを比較してみます。
シクロヘキサンの立体配座 より
初期座標の構造用意
いす形はPubChemで公開されている構造を使用します。
3D structure image of CID 8078 Cyclohexane (Compound)
PubChem Identifier: CID 8078 URL: https://pubchem.ncbi.nlm.nih.gov/compound/8078#section=3D-Conformer
[7]:
from ase.io import read
chair_atoms = read("../input/Conformer3D_CID_8078.sdf")
view_ngl(chair_atoms, representations=["ball+stick"], w=400, h=300)
舟形は今回は、いす形から座標値を手動で変更して生成してみましょう。
※: コードの意味はわからなくてもスキップしていただいて構いません。(原子を180度回転させたものをつくり、原子2, 10, 11のみ回転させた座標のものに上書きしています。)
[8]:
boat_atoms = chair_atoms.copy()
axis = boat_atoms.positions[[1, 5]].mean(axis=0) - boat_atoms.positions[[0, 4]].mean(axis=0)
boat_atoms.rotate(180, axis)
# inds = [1, 3, 5, 8, 9, 12, 13, 16, 17]
inds = [2, 10, 11]
boat_atoms.positions[inds] = chair_atoms.positions[inds]
[9]:
view_ngl(boat_atoms, representations=["ball+stick"], w=400, h=300)
[10]:
chair_atoms.calc = calculator
opt = LBFGS(chair_atoms, trajectory="output/lbfgs_cyclohexane_chair.traj")
opt.run()
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 02:20:03 -76.855956* 0.1778
LBFGS: 1 02:20:03 -76.858225* 0.1151
LBFGS: 2 02:20:03 -76.860049* 0.0998
LBFGS: 3 02:20:03 -76.862485* 0.0499
[10]:
True
[11]:
boat_atoms.calc = calculator
opt = LBFGS(boat_atoms, trajectory="output/lbfgs_cyclohexane_boat.traj")
opt.run()
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 02:20:14 -71.025783* 5.3644
LBFGS: 1 02:20:14 -72.330446* 4.4228
LBFGS: 2 02:20:14 -74.196692* 3.2822
LBFGS: 3 02:20:14 -75.112910* 4.2950
LBFGS: 4 02:20:14 -75.580790* 2.7280
LBFGS: 5 02:20:14 -75.743022* 2.4649
LBFGS: 6 02:20:15 -75.854447* 1.8966
LBFGS: 7 02:20:15 -76.089088* 0.7839
LBFGS: 8 02:20:15 -76.161671* 0.6727
LBFGS: 9 02:20:15 -76.249390* 0.8271
LBFGS: 10 02:20:15 -76.300558* 0.7451
LBFGS: 11 02:20:15 -76.366199* 0.6365
LBFGS: 12 02:20:15 -76.410962* 0.5463
LBFGS: 13 02:20:15 -76.432972* 0.2407
LBFGS: 14 02:20:15 -76.442167* 0.2395
LBFGS: 15 02:20:15 -76.450881* 0.2134
LBFGS: 16 02:20:15 -76.458265* 0.1773
LBFGS: 17 02:20:15 -76.462880* 0.1930
LBFGS: 18 02:20:15 -76.466240* 0.1989
LBFGS: 19 02:20:15 -76.470400* 0.2080
LBFGS: 20 02:20:15 -76.476647* 0.2489
LBFGS: 21 02:20:15 -76.483845* 0.2404
LBFGS: 22 02:20:15 -76.489993* 0.1807
LBFGS: 23 02:20:15 -76.494768* 0.2158
LBFGS: 24 02:20:15 -76.499497* 0.2090
LBFGS: 25 02:20:16 -76.504614* 0.1982
LBFGS: 26 02:20:16 -76.509015* 0.2078
LBFGS: 27 02:20:16 -76.512713* 0.2074
LBFGS: 28 02:20:16 -76.517285* 0.1909
LBFGS: 29 02:20:16 -76.523674* 0.2553
LBFGS: 30 02:20:16 -76.530986* 0.2316
LBFGS: 31 02:20:16 -76.537306* 0.2286
LBFGS: 32 02:20:16 -76.542625* 0.2381
LBFGS: 33 02:20:16 -76.548161* 0.2532
LBFGS: 34 02:20:16 -76.554129* 0.2480
LBFGS: 35 02:20:16 -76.558687* 0.1736
LBFGS: 36 02:20:16 -76.560842* 0.0990
LBFGS: 37 02:20:16 -76.561720* 0.0792
LBFGS: 38 02:20:16 -76.562696* 0.0877
LBFGS: 39 02:20:16 -76.563681* 0.0767
LBFGS: 40 02:20:16 -76.564258* 0.0412
[11]:
True
構造最適化後の構造を可視化して確認してみましょう。
[12]:
view_ngl([chair_atoms, boat_atoms], representations=["ball+stick"], w=400, h=300, replace_structure=True)
また、それぞれのエネルギー差を比較してみます。
エネルギーの単位はeVですが、単位の変換する場合、units module を用いることができます。
[13]:
from ase.units import Bohr,Rydberg,kJ,kB,fs,Hartree,mol,kcal
E_chair = chair_atoms.get_potential_energy()
E_boat = boat_atoms.get_potential_energy()
E_diff = E_boat - E_chair
print(f"E_chair {E_chair:.2f} eV")
print(f"E_boat {E_boat:.2f} eV")
print(f"E_diff {E_diff:.2f} eV = {E_diff / (kJ/mol):.2f} kJ/mol")
E_chair -76.86 eV
E_boat -76.56 eV
E_diff 0.30 eV = 28.77 kJ/mol
いす形の方が舟形よりもエネルギーが低く、安定であることがわかりました。 また、そのエネルギー差も上図で示されていた25kJ/mol に近い値を得ることができました。
このように構造最適化では、局所安定点を求めるため、初期構造が重要となってきます。 (おなじシクロヘキサンでも初期構造の用意の仕方によって構造最適化後の結果が異なります。)
Opt Crystals¶
次に周期構造のある系に対する構造最適化を行ってみます。 上記の周期構造がない有機分子では各原子の座標のみが最適化対象でしたが、周期系ではセルサイズ も最適化対象に加わります。
また、周期構造を持つ結晶系 (Bulk)では、その結晶構造によりエネルギーが異なります。
1. Si結晶¶
ここではSiの結晶構造として、どのような構造がエネルギー安定となるのかを調べてみましょう。 結晶構造の候補として、以下のような構造を試してみます。
sc: 単純立方格子 (Simple Cubic)
bcc: 体心立方格子 (Body Centered Cubic)
fcc: 面心立方格子 (Face Centered Cubic)
hcp: 六方最密充填構造 （Hexagonal Close Packed)
diamond: ダイヤモンド構造
エネルギーは１原子あたりのエネルギーを算出します。
まず、Si のBCC構造に対し、その格子定数 a を様々な値に変えた場合のポテンシャルエネルギー及びその時の密度 (density) を計算するコードは以下のようになります。
ここでは、np.linspace 関数を用いて格子定数 a を 2.5A から7.0A の間で5つ生成し、それぞれの格子定数におけるエネルギーを求めています。
[14]:
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL_U0, model_version="v3.0.0")
calculator = ASECalculator(estimator)
[15]:
from ase.build import bulk
a_array = np.linspace(2.5, 7.0, 5)
bcc_energy_list = []
bcc_density_list = []
for a in a_array:
atoms = bulk("Si", "bcc", a=a)
atoms.calc = calculator
# calculate potential energy for 1 atom (eV/atom)
E_pot = atoms.get_potential_energy() / len(atoms)
density = len(atoms) / atoms.get_volume()
bcc_energy_list.append(E_pot)
bcc_density_list.append(atoms.get_volume())
print("lattice a: ", list(a_array))
print("Density : ", bcc_density_list)
print("Energy : ", bcc_energy_list)
lattice a: [2.5, 3.625, 4.75, 5.875, 7.0]
Density : [7.8125, 23.817382812499996, 53.58593750000001, 101.38964843749996, 171.50000000000006]
Energy : [-1.1319000978713731, -3.3385418260066886, -1.325669089403413, -0.1946264644375943, -0.007487893090558714]
今回は、様々な構造に対して計算を行うため、これを関数化してみましょう。
[16]:
from ase.build import bulk
def calc_density_energy(calculator, symbol: str = "Fe", crystalstructure: str = "bcc"):
a_array = np.linspace(2.0, 8.0, 30)
energy_list = []
density_list = []
for a in a_array:
atoms = bulk(symbol, crystalstructure, a=a)
atoms.calc = calculator
E_pot = atoms.get_potential_energy() / len(atoms)
density = len(atoms) / atoms.get_volume()
energy_list.append(E_pot)
density_list.append(density)
return np.array(energy_list), np.array(density_list)
この calc_density_energy 関数を用いることで、様々な構造に対するエネルギーと密度のリストを得ることができます。
[17]:
bcc_energy, bcc_density = calc_density_energy(calculator, "Si", "bcc")
fcc_energy, fcc_density = calc_density_energy(calculator, "Si", "fcc")
sc_energy, sc_density = calc_density_energy(calculator, "Si", "sc")
hcp_energy, hcp_density = calc_density_energy(calculator, "Si", "hcp")
diamond_energy, diamond_density = calc_density_energy(calculator, "Si", "diamond")
計算されたエネルギー値を見てみると、aの値が小さすぎる場合 (現実的でないほど密度を大きくしてしまった場合)にはエネルギーの値が不安定に出力されてしまっていることがわかります。
[18]:
bcc_energy
[18]:
array([-1.51431809e+03, 1.90312736e+02, 1.00883824e+00, -2.69269530e+00,
-3.69642360e+00, -4.01363424e+00, -3.93359487e+00, -3.67273086e+00,
-3.28083585e+00, -2.88532212e+00, -2.49506160e+00, -2.09706308e+00,
-1.72758801e+00, -1.41201213e+00, -1.11850814e+00, -8.32015884e-01,
-5.81543805e-01, -3.93798904e-01, -2.63483682e-01, -1.65962615e-01,
-9.01755903e-02, -5.25514094e-02, -4.20469155e-02, -3.06421093e-02,
-7.48789309e-03, -7.48789309e-03, -7.48789309e-03, -7.48789309e-03,
-7.48789309e-03, -7.48789309e-03])
以下はこのような値は避けて、プロットを行っています。
np.where関数を使って、絶対値が5以下になっているエネルギーをとっているindex をinds としてとりだし、a_array[inds]でハズレ値を除いたプロットを行っています。
[19]:
a_array = np.linspace(2.0, 8.0, 30)
inds = np.where(np.abs(sc_energy) < 5.0)[0]
plt.plot(a_array[inds], sc_energy[inds], label="sc")
inds = np.where(np.abs(bcc_energy) < 5.0)[0]
plt.plot(a_array[inds], bcc_energy[inds], label="bcc")
inds = np.where(np.abs(fcc_energy) < 5.0)[0]
plt.plot(a_array[inds], fcc_energy[inds], label="fcc")
inds = np.where(np.abs(hcp_energy) < 5.0)[0]
plt.plot(a_array[inds], hcp_energy[inds], label="hcp")
inds = np.where(np.abs(diamond_energy) < 5.0)[0]
plt.plot(a_array[inds], diamond_energy[inds], label="diamond")
plt.xlabel("a [A]")
plt.ylabel("Energy per atom [eV]")
plt.legend()
plt.title("Si bulk energy for various structure")
plt.ylim(-6.0, 0)
#plt.yscale("log")
plt.show()
上図は格子定数aを横軸にとっていますが、この場合は構造ごとにエネルギーが小さくなる位置が異なり比較しづらい事がわかります。
下図では、密度を横軸にとって、エネルギーをプロットしてみます。
[20]:
inds = np.where(np.abs(sc_energy) < 5.0)[0]
plt.plot(sc_density[inds], sc_energy[inds], label="sc")
inds = np.where(np.abs(bcc_energy) < 5.0)[0]
plt.plot(bcc_density[inds], bcc_energy[inds], label="bcc")
inds = np.where(np.abs(fcc_energy) < 5.0)[0]
plt.plot(fcc_density[inds], fcc_energy[inds], label="fcc")
inds = np.where(np.abs(hcp_energy) < 5.0)[0]
plt.plot(hcp_density[inds], hcp_energy[inds], label="hcp")
inds = np.where(np.abs(diamond_energy) < 5.0)[0]
plt.plot(diamond_density[inds], diamond_energy[inds], label="diamond")
plt.xlabel("Density [1/A^3]")
plt.ylabel("Energy per atom [eV]")
plt.legend()
plt.title("Si bulk energy for various structure")
plt.ylim(-6.0, 0)
#plt.yscale("log")
plt.show()
密度を横軸に取ることで見やすい図になりました。
このように、格子定数を変えながら計算をすることで、SC, BCC, FCC, HCP, Diamondそれぞれの構造でエネルギーが最小になりそうな格子定数を求めることができます。 今回の計算結果で、それぞれの構造に対する最低エネルギーを比較することで、Si はDiamondが安定であるということが分析できました。
2. Fe結晶の格子定数最適化¶
次は、この格子定数を構造最適化アルゴリズムを用いることにより、効率よく求めてみましょう。
周期構造のある系の最適化を行う際には、原子座標のみを最適化するのか、Cell sizeも含めて最適化するのかによりその方法が異なります。 Cell も最適化の対象としたい場合、ASE のatomsに対して、以下に示すような Filter クラスを適用する必要があります。
Table: 最適化対象と、使用するFilter
Position |
Cell |
Filter |
|---|---|---|
✓ |
不要 (Default挙動) |
|
✓ |
|
|
✓ |
✓ |
|
座標は固定して、CellのみをOpt –>
StrainFilterCell固定して+座標のみをOpt –> Filter未使用
Cell+座標両方をOpt –>
ExpCellFilter
まずはじめにStrainFilterを使う例を紹介します。
ここではFe原子のBCC構造・FCC構造の格子定数を最適化してみます。 BCC構造・FCC構造はprimitive cellに原子が１つしか含まれないため、原子座標を最適化する必要はなく、格子定数のみの最適化を行えば十分なため、StrainFilterを用いています。
[21]:
from ase.constraints import ExpCellFilter, StrainFilter
fe_bcc_atoms = bulk("Fe", "bcc", a=3.0)
print("Before opt:", fe_bcc_atoms.cell)
fe_bcc_atoms.calc = calculator
fe_bcc_atoms = StrainFilter(fe_bcc_atoms)
opt = LBFGS(fe_bcc_atoms, trajectory="output/lbfgs_fe_bcc.traj")
opt.run()
Before opt: Cell([[-1.5, 1.5, 1.5], [1.5, -1.5, 1.5], [1.5, 1.5, -1.5]])
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 02:20:35 -4.865733* 2.2558
LBFGS: 1 02:20:35 -4.933622* 2.0702
LBFGS: 2 02:20:35 -4.434309* 8.0541
LBFGS: 3 02:20:35 -4.992679* 0.5748
LBFGS: 4 02:20:35 -4.996655* 0.1459
LBFGS: 5 02:20:35 -4.996941* 0.0051
[21]:
True
[22]:
from ase.constraints import ExpCellFilter, StrainFilter
fe_fcc_atoms = bulk("Fe", "fcc", a=3.0)
print("Before opt:", fe_fcc_atoms.cell)
fe_fcc_atoms.calc = calculator
fe_fcc_atoms = StrainFilter(fe_fcc_atoms)
opt = LBFGS(fe_fcc_atoms, trajectory="output/lbfgs_si_fcc.traj")
opt.run()
Before opt: Cell([[0.0, 1.5, 1.5], [1.5, 0.0, 1.5], [1.5, 1.5, 0.0]])
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 02:20:36 -2.789322* 16.4962
LBFGS: 1 02:20:36 -4.806755* 3.0130
LBFGS: 2 02:20:36 -4.878337* 0.6926
LBFGS: 3 02:20:36 -4.882889* 0.1071
LBFGS: 4 02:20:36 -4.883017* 0.0137
[22]:
True
[23]:
fe_bcc_atoms.atoms.cell
[23]:
Cell([[-1.4243732594602483, 1.4243732644003524, 1.4243732505018425], [1.4243732675507748, -1.4243732687974981, 1.4243732629895147], [1.42437325401479, 1.4243732633520398, -1.4243732548989883]])
[24]:
fe_fcc_atoms.atoms.cell
[24]:
Cell([[-2.272649880573051e-09, 1.7255768493652164, 1.7255767924285237], [1.7255767836800167, 1.862583482829644e-09, 1.7255767882932904], [1.7255767714053183, 1.7255768329552847, 1.4137281709826914e-08]])
StrainFilterを用いて格子定数の最適化を行った場合、結晶のa, b, c軸の対称性を保つ考慮がされないためそれぞれ値が少し違った値となってしまっています。 実際は後述するExpCellFilter や、FixSymmetry を活用することで、対称性を保ったままの最適化も可能です。
ここでは、最終的な格子定数を計算するために平均値をとることにします。 ASEのbulk methodに指定する引数 a に対して、BCC構造のセルは
[[-a/2, a/2, a/2],
[a/2, -a/2, a/2],
[a/2, a/2, -a/2]]
FCC構造のセルは
[[0, a/2, a/2],
[a/2, 0, a/2],
[a/2, a/2, 0]]
となっていることを用いて、下記の計算を行っています。
[25]:
a = np.abs(fe_bcc_atoms.atoms.cell).mean() * 2
E_pot = fe_bcc_atoms.get_potential_energy()
print(f"BCC opt: a = {a:.2f} A, E = {E_pot:.2f} eV")
a = fe_fcc_atoms.atoms.cell.sum() / 3
E_pot = fe_fcc_atoms.get_potential_energy()
print(f"FCC opt: a = {a:.2f} A, E = {E_pot:.2f} eV")
BCC opt: a = 2.85 A, E = -5.00 eV
FCC opt: a = 3.45 A, E = -4.88 eV
構造最適化の結果として、BCC構造では a=2.85 A のときにエネルギーが最小となる事がわかりました。
また、FeはBCC構造のほうがFCC構造よりも安定であるという、実験結果と整合する結果を得ることができました。
構造最適化の途中経過を可視化してみると、原子座標は変わらずに、格子定数のみが最適化されていることが確認できます。
[26]:
view_ngl(Trajectory("output/lbfgs_fe_bcc.traj"), w=400, h=300, replace_structure=True)
StrainFilterを使う場合の注意点として、セルのa, b, c軸が独立して最適化されてしまうため、その比率を保ちたい場合にも保たれない可能性があります。
3. Si結晶diamond 構造の格子定数最適化¶
次に、Si結晶のdiamond構造に対して、最適化を行ってみます。 ここでは、原子座標及びセルサイズ双方を最適化するExpCellFilter を使用してみます。
ExpCellFilterでは、hydrostatic_strain=TrueというOptionをつけることによって、セルのa, b, c軸長さの比率を保ったまま最適化が行なえます。
[27]:
from ase.constraints import ExpCellFilter, StrainFilter
si_diamond_atoms = bulk("Si", "diamond")
si_diamond_atoms.calc = calculator
print("cell size before: ", si_diamond_atoms.cell)
si_diamond_atoms = ExpCellFilter(si_diamond_atoms, hydrostatic_strain=True)
opt = LBFGS(si_diamond_atoms, trajectory="output/lbfgs_si_diamond.traj")
opt.run()
print("cell size after : ", si_diamond_atoms.atoms.cell)
cell size before: Cell([[0.0, 2.715, 2.715], [2.715, 0.0, 2.715], [2.715, 2.715, 0.0]])
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 02:20:49 -9.108032* 0.4967
LBFGS: 1 02:20:49 -9.113262* 0.0002
cell size after : Cell([[0.0, 2.734332539572064, 2.734332539572064], [2.734332539572064, 0.0, 2.734332539572064], [2.734332539572064, 2.734332539572064, 0.0]])
cell sizeのbefore, after を見比べると、セルサイズの比率が保たれたまま、構造緩和を行えたことがわかります。
[28]:
view_ngl(Trajectory("output/lbfgs_si_diamond.traj"), w=400, h=300, replace_structure=True)
4. SiO\(_2\) 結晶¶
これまでは1元素のみを含む結晶を扱ってみました。最後に実用的な例として、SiO\(_2\)結晶についての最適化を行ってみます。
SiO\(_2\)は、現実でもQuartz (石英), Cristobalite, Tridymite (鱗珪石)といった、様々な構造を取ることが知られています。 このように同じ組成で多くの結晶系が存在することを多形と呼びます。
今回は、Quartz構造をMaterials projectから用意し、構造最適化を行っています。 座標値・セルサイズ双方を同時に最適化するため、ExpCellFilter を用いています。
[29]:
atoms = read("../input/SiO2_mp-6930_conventional_standard.cif")
atoms.calc = calculator
view_ngl(atoms, representations=["ball+stick"])
[30]:
opt = LBFGS(ExpCellFilter(atoms), trajectory="output/lbfgs_sio2.traj")
opt.run()
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 02:20:57 -57.136945* 0.0379
[30]:
True
[31]:
view_ngl(Trajectory("output/lbfgs_sio2.traj"), representations=["ball+stick"])
[Note] 構造最適化はCalculator (Potential energyの関数系)が変わる度に必要
Materials projectの構造も、構造最適化を行った構造であるため、すでに完全な安定構造かと思うかもしれません。 しかし、厳密な局所点(極小のエネルギーを取る構造)は、ポテンシャルエネルギーが変わると異なります。
[コラム] 構造最適化計算と現実世界¶
この章で紹介した構造最適化計算は、今後紹介するほぼすべての原子シミュレーションにとって土台となる計算となります。しかし、構造最適化計算の内部で使われている計算手法はいくらか技巧的なものであり、このチュートリアルではじめて原子シミュレーションに触れた人にとっては直感的には必要性がわかりにくいかもしれません。そこでこのコラムでは、現実世界との対応関係という観点から構造最適化計算を見直してみます。
構造最適化の必要性¶
現実世界の物質の原子配置がどのように決まっているかについて考えてみますと、一つ一つの原子は特に拘束もなく自由に動き回っているので、余計なエネルギーは散逸し、エネルギーの低い安定な構造へと変化していきます。一旦温度の影響を無視して十分に温度が低い場合を考えてみると、その場合は原子の運動エネルギーは0になり、ポテンシャルエネルギーは極小値をとることが考えられます。その意味では構造最適化計算は現実世界における絶対零度の状態を考えていることに対応しています。
有限温度の場合にはどうでしょうか。統計力学の理論から、平衡状態にある系では温度と運動エネルギーは一定の関係にあることが知られており、大雑把な計算では1 Kの上昇はおおまかに1/10000 eV/atomの運動エネルギーの上昇に相当します。逆にいうと、1 eV/atomのエネルギーの上昇は温度に換算して10000 K程度の構造に相当するということです。このように有限温度の場合でもしばしばエネルギーを最適化した構造が十分に良い近似となることが多いため、シミュレーションのはじめにまずは構造最適化計算をしよう(有限温度の場合にはその後で温度を上昇させよう)という手順をとることが一般に行われています。
構造最適化計算を行わずにシミュレーションをしてしまうと、前提条件の違いにより非現実的な結果をもたらしてしまう可能性があります。例えば鉄の結晶構造を1%小さくとった場合、原子構造を可視化しただけであれば特段の変化は見えないかもしれません。しかし現実世界で鉄を1%縮める操作を考えてみればわかるように、これは非常に強い圧縮力(具体的には1 GPaを超える圧力)のかかった鉄に相当するので、シミュレーションの結果を変えてしまうことが想定されます。
構造最適化手法と現実世界の構造緩和¶
構造最適化計算の種類がたくさんあることについて、現実世界との対応から考えてみます。適当な温度の高い材料、例えば溶けた鉄について、ここから熱を取り除いて少し冷えていくときに起きる現象を考えてみましょう。鉄原子は熱によって常に振動しており、様々な原子配置をとりながら、少しずつポテンシャルエネルギーの低い構造へと遷移していきます。最終的にはエネルギーの低い、原子がきれいに整列した結晶構造へと遷移していきます。この意味では現実世界における構造最適化は焼きなまし法(SA: Simulated Annealing)と呼ばれる最適化手法の一種とみなすことができます。(歴史的には、こちらが焼きなまし法のアイデア元に相当します。)
これは原子が長時間運動した結果として得られるものです。つまり、原理的には原子シミュレーションでも、十分に長い分子動力学計算が行えれば現実世界と同じように安定な構造が得られることを意味しています。
しかし、ここで問題になるのが時間スケールのギャップです。原子の動きは典型的にはフェムト秒(\(10^{-15}\) s)オーダーの現象ですが、現実世界の材料を考えるときは一般にはるかに長い時間スケールが登場します。このため、現実世界の構造を原子スケールシミュレーションの言葉で書き直すと、\(10^{15}\)やそれ以上のスケールでの極めて長い時間をかけた構造最適化計算の結果とみなすことができます。また、焼きなまし法はロバスト性の高い最適化手法であり、十分に長い時間をかけることでより大域的な最適解へとたどり着くことが知られています。このため、現実に存在する物質は実質的に大域最適化問題が解かれた原子構造が得られていると考えることができます。
このようなことから、原子シミュレーションでも短時間でより大域的な最適解を探すことができる手法が開発されてきました。これは一般に計算時間とのトレードオフとなるため、ケースによって使い分けられます。結晶構造のように目的の構造が既にわかっていて十分に近い構造を初期構造としてとれる場合には、上の例にあるBFGS法のように局所最適解を求める手法が好まれます。逆にアモルファス構造のようによくわからない構造を対象とする場合には、より大域的な最適解を求めるためにより注意が向けられます。(より現実的になるよう、ゆっくりと温度を下げていくメルトクエンチ法などが使われます。)ASEを含む原子シミュレーションのライブラリに様々な最適解手法が搭載されているのはこのためです。
余談ですが、実は現実世界でも大域最適解になりきらないケースは存在します。例えば炭素の高圧相であるダイヤモンドは地上では安定ではなく、黒鉛(グラファイト)のほうがより安定なのですが、構造変化に必要なエネルギーが高いため構造が変化する速度が極めて遅く、人間の時間スケールでは実質的に安定な構造(準安定状態)として扱われます。製造時に積極的に利用している有名な例としては鉄の焼入れがあります。日本刀を作るときに熱した鉄を水に漬けることで硬くする技術が有名ですが、温度を急冷させることで構造緩和にかかる時間を短縮し、安定構造に遷移することを防いでいます。これによってより望ましい強度を持つ準安定状態を得ています。 このように、相変態や化学変化を伴う現象には日常スケールでも緩慢に起きる場合が多いです。これは原子シミュレーションでも同様で、相変態や化学反応を見たい場合には単純に長時間の動力学計算を行うよりも効率的な計算手法が使われることがあります。ひとつの具体例として、5章で説明があるNEB法があります。