Appendix 2: Opt algorithm fmax experiment

本節では構造最適化のパラメーターであるfmaxについて考察します。

!pip install more_itertools

from ase import Atoms
from ase.constraints import ExpCellFilter, FixAtoms
from ase.build import bulk, molecule, add_adsorbate, surface
from ase.optimize import FIRE
from ase.units import kB
import pandas as pd
from tqdm.auto import tqdm
import numpy as np
from more_itertools import windowed
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
from pfcc_extras.structure.ase_rdkit_converter import smiles_to_atoms, atoms_to_smiles
from pfcc_extras.visualize.view import view_ngl
from matlantis_features.features.common.fire_lbfgs import FIRELBFGS

/home/jovyan/.local/lib/python3.7/site-packages/matlantis_features/features/common/fire_lbfgs.py:12: FutureWarning: matlantis_features.features.common.fire_lbfgs is deprecated; use ase_ext.optimize
  FutureWarning,

def max_distance(atoms1: Atoms, atoms2: Atoms) -> float:
    return float(np.max(np.linalg.norm(atoms1.positions - atoms2.positions, axis=1)))

calc_mol = ASECalculator(Estimator(calc_mode=EstimatorCalcMode.MOLECULE, model_version="v2.0.0"))
calc = ASECalculator(Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v2.0.0"))
calc_d3 = ASECalculator(Estimator(calc_mode=EstimatorCalcMode.CRYSTAL_PLUS_D3, model_version="v2.0.0"))

fmax=0.05eV/Aで十分な例

ASEにおいてfmaxのデフォルト値は0.05eV/Aですが、これは多くの場合において十分な精度を持っています。

分子と固体の例でそれぞれ確認してみましょう。

分子の例

fmax_list = [0.05, 0.01, 0.001]

atoms = molecule("C2H6")
atoms.rattle(stdev=0.1)
atoms.calc = calc_mol
images = []
for fmax in tqdm(fmax_list):
    with FIRE(atoms, logfile=None) as opt:
        opt.run(fmax=fmax)
        images.append(atoms.copy())

for atoms in images:
    atoms.calc = calc_mol

results = {
    "fmax": [],
    "energy (eV)": [],
    "distance (A)": [],
}
for (fmax_old, old), (fmax_new, new) in windowed(zip(fmax_list, images), 2):
    de = new.get_potential_energy() - old.get_potential_energy()
    dx = max_distance(new, old)
    results["energy (eV)"].append(de)
    results["distance (A)"].append(dx)
    results["fmax"].append(f"fmax {fmax_old}->{fmax_new}")
pd.DataFrame(results)

	fmax	energy (eV)	distance (A)
0	fmax 0.05->0.01	-0.002484	0.030840
1	fmax 0.01->0.001	-0.000402	0.017346

固体の例

fmax_list = [0.05, 0.01, 0.001]

atoms = bulk("Pt") * (4, 4, 4)
atoms.rattle(stdev=0.1)
atoms.calc = calc_mol
images = []
for fmax in tqdm(fmax_list):
    with FIRE(atoms, logfile=None) as opt:
        opt.run(fmax=fmax)
        images.append(atoms.copy())

for atoms in images:
    atoms.calc = calc_mol

results = {
    "fmax": [],
    "energy (eV)": [],
    "distance (A)": [],
}
for (fmax_old, old), (fmax_new, new) in windowed(zip(fmax_list, images), 2):
    de = new.get_potential_energy() - old.get_potential_energy()
    dx = max_distance(new, old)
    results["energy (eV)"].append(de)
    results["distance (A)"].append(dx)
    results["fmax"].append(f"fmax {fmax_old}->{fmax_new}")
pd.DataFrame(results)

	fmax	energy (eV)	distance (A)
0	fmax 0.05->0.01	-0.005323	0.016975
1	fmax 0.01->0.001	-0.001351	0.013533

いずれの場合も、最も動いた原子でも0.01A程度しか動かず、エネルギーも0.01eV程度しか変わりません。

0.01eVというと、300Kで存在比が1.5倍になる程度のエネルギー差です。常温で扱う材料を探索したい場合、この値は数値計算で欲しい精度を超えています。

また、実験化学の精度に対応するとされるchemical accuracyは0.04eV程度ですが、多くの数値計算手法はその精度に達していないのが現状であり、その意味でもこの精度は必要ないことが多いでしょう。

0.01eVのエネルギー差がある2つの構造の300Kにおける存在比はボルツマン因子から以下のように求めます。

np.exp(0.01 / (300 * kB))

1.4722876299701386

fmax=0.05eVで不十分な例

エネルギーに関して精度不足な例

少し意地悪な例ですが、2,3-Dimethyl-2-buteneがPt表面にゆるやかに吸着しているような例を考えてみましょう。

このような例ではゆるやかな力が働くため、fmax=0.05eV/Aでは構造最適化が途中で止まってしまい、さらに最適化した場合と比べて大きなエネルギー差が生じてしまいます。

このような状況が発生する可能性があるならば、なるべく小さなfmaxを使ったほうがより信頼のおける構造最適化が出来るでしょう。

bulk1111 = bulk("Pt")
bulk1111.calc = calc_d3
with FIRELBFGS(ExpCellFilter(bulk1111), logfile=None) as opt:
    opt.run(0.0001)

atoms = surface(bulk1111, (1, 1, 1), 4, vacuum=20.0) * (4, 4, 1)
c = atoms.cell[2, 2] / 2
atoms.constraints = [FixAtoms(mask=[atom.position[2] < c for atom in atoms])]
atoms.rattle(stdev=0.1)
ads = smiles_to_atoms("CC(=C(C)C)C")
v = (ads.positions[0] - ads.positions[1]) / 2
p = (atoms.cell[0, :2] + atoms.cell[1, :2]) / 2
add_adsorbate(atoms, ads, 4.0, position=tuple(p))
atoms.positions[64:] = atoms.positions[64:] + v
atoms.calc = calc_d3

view_ngl(atoms, representations=["ball+stick"])

fmax_list = [0.05, 0.01, 0.001]

images = []
for fmax in tqdm(fmax_list):
    with FIRELBFGS(atoms, logfile=None) as opt:
        opt.run(fmax=fmax)
        images.append(atoms.copy())

for atoms in images:
    atoms.calc = calc_mol

results = {
    "fmax": [],
    "energy (eV)": [],
    "distance (A)": [],
}
for (fmax_old, old), (fmax_new, new) in windowed(zip(fmax_list, images), 2):
    de = new.get_potential_energy() - old.get_potential_energy()
    dx = max_distance(new, old)
    results["energy (eV)"].append(de)
    results["distance (A)"].append(dx)
    results["fmax"].append(f"fmax {fmax_old}->{fmax_new}")
pd.DataFrame(results)

	fmax	energy (eV)	distance (A)
0	fmax 0.05->0.01	-0.213606	4.068967
1	fmax 0.01->0.001	-0.987591	2.578896

view_ngl(images, representations=["ball+stick"], replace_structure=True)

座標に関して精度不足な例（１）

Pt(111)面にエチレンが吸着したような例では、fmax=0.05eVでは中途半端な構造で構造最適化が止まってしまいます。

bulk1111 = bulk("Pt")
bulk1111.calc = calc
with FIRELBFGS(ExpCellFilter(bulk1111), logfile=None) as opt:
    opt.run(0.0001)

atoms = surface(bulk1111, (1, 1, 1), 4, vacuum=20.0) * (4, 4, 1)
c = atoms.cell[2, 2] / 2
atoms.constraints = [FixAtoms(mask=[atom.position[2] < c for atom in atoms])]
atoms.calc = calc
atoms.rattle(stdev=0.1)
add_adsorbate(atoms, smiles_to_atoms("C=C"), 2.0, position=tuple((atoms.cell[0, :2] + atoms.cell[1, :2]) / 2))

view_ngl(atoms, representations=["ball+stick"])

images = []
for fmax in tqdm(fmax_list):
    with FIRELBFGS(atoms, logfile=None) as opt:
        opt.run(fmax=fmax)
        images.append(atoms.copy())

for atoms in images:
    atoms.calc = calc_mol

results = {
    "fmax": [],
    "energy (eV)": [],
    "distance (A)": [],
}
for (fmax_old, old), (fmax_new, new) in windowed(zip(fmax_list, images), 2):
    de = new.get_potential_energy() - old.get_potential_energy()
    dx = max_distance(new, old)
    results["energy (eV)"].append(de)
    results["distance (A)"].append(dx)
    results["fmax"].append(f"fmax {fmax_old}->{fmax_new}")
pd.DataFrame(results)

	fmax	energy (eV)	distance (A)
0	fmax 0.05->0.01	-0.020390	0.285722
1	fmax 0.01->0.001	-0.001203	0.083779

view_ngl(images, representations=["ball+stick"], replace_structure=True)

中途半端で止まってしまった構造と最適化しきった構造を見比べると違った印象をうけると思います。

とはいえこの例ではエネルギー差は大きくないため、場合によってはこの違いを無視するという判断を行っても良いかもしれません。

座標について精度不足な例（２）

構造最適化が難しい系として良く知られているトルエンのメチル基の回転を試してみます。

atoms = smiles_to_atoms("Cc1ccccc1")
tmp = atoms[7:10]
tmp.rotate([1.0, 0.0, 0.0], 15.0)
atoms.positions[7:10] = tmp.positions
atoms.calc = calc_mol

view_ngl(atoms, representations=["ball+stick"])

images = []
for fmax in tqdm(fmax_list):
    with FIRELBFGS(atoms, logfile=None) as opt:
        opt.run(fmax=fmax)
        images.append(atoms.copy())

for atoms in images:
    atoms.calc = calc_mol

results = {
    "fmax": [],
    "energy (eV)": [],
    "distance (A)": [],
}
for (fmax_old, old), (fmax_new, new) in windowed(zip(fmax_list, images), 2):
    de = new.get_potential_energy() - old.get_potential_energy()
    dx = max_distance(new, old)
    results["energy (eV)"].append(de)
    results["distance (A)"].append(dx)
    results["fmax"].append(f"fmax {fmax_old}->{fmax_new}")
pd.DataFrame(results)

	fmax	energy (eV)	distance (A)
0	fmax 0.05->0.01	-0.001788	0.028956
1	fmax 0.01->0.001	-0.001723	0.260170

view_ngl(images, representations=["ball+stick"], replace_structure=True)

トルエンのメチル基の回転は角度が変わってもエネルギー差があまり変わらないことが判ると思います。

このような例では小さなfmaxを用いないと最安定構造の予測が出来ません。 0.002eVのエネルギー差は300Kで存在比が1.08倍になる程度の影響しかないので、ほとんど影響がないと言っても構わないでしょう。

従ってこのような小さなエネルギー差は場合によっては無視してかまいません。

ただし、振動解析のようなエネルギーが極小であることを前提とした解析を行いたい場合は時として解析に悪影響を及ぼすかもしれません。

np.exp(0.002 / (300 * kB))

1.080434721876578