Optimization algorithm¶
前節まででは様々な事例に対して、 構造最適化を適用してみました。 本節では、構造最適化の際に適用を行った局所最適化アルゴリズムについて学んでいきます。
アルゴリズムの種類¶
局所最適化アルゴリズムとして、ASEではMolecular Dynamics (MD)に似た振る舞いで最適化を行う FIRE や、2次関数近似をおこなって最適化を行う BFGS, LBFGS, BFGSLineSearch=QuasiNewton, LBFGSLineSearch法などが提供されています。
アルゴリズム |
グループ |
特徴 |
|---|---|---|
MDMin |
MD-like |
MDを行うが、運動量とForceの内積が負になると運動量を0にリセットする |
FIRE |
MD-like |
MDを行うが、様々な追加の工夫を入れて高速かつ頑健にしたもの |
BFGS |
準Newton法 |
最適化のトラジェクトリからヘシアンを近似し、近似ヘシアンを用いてNewton法のアルゴリズムを実行 |
LBFGS |
準Newton法 |
BFGS法を高速・低メモリで動くようにしたもの |
BFGSLineSearch |
準Newton法 |
BFGSで最適化のステップの方向を決定し、ステップの大きさはLineSearchで決定 |
LBFGSLineSearch |
準Newton法 |
LBFGSで最適化のステップの方向を決定し、ステップの大きさはLineSearchで決定 |
MDMinはハイパーパラメーター依存性が大きく、構造最適化に失敗することも多いためベンチマークから除外し、他のOptimizerを見ていきます。
FIRE¶
FIREでは基本的にMDを行いますが、高速に収束するように様々な工夫が加えられています。
ニュートン法¶
ニュートン法ではヘシアンと勾配を用いて次のステップを決定します。ニュートン法は二次収束する手法であり、ポテンシャルが二次関数に近い時（極小点に近い時）真の解に二乗の速度で近づいていき非常に高速です。 ヘシアンを \(H\)、勾配を \(\vec{g}\) として、次の最適化ステップ \(\vec{p}\) を以下の式で表現します。
\(\vec{p} = - H^{-1} \vec{g}\)
仮にポテンシャルが厳密に二次関数であれば1ステップで極小値に収束します。
準ニュートン法¶
ヘシアンの評価は非常に重い（例えば数値計算なら6N回Forceを計算して求めることが出来る）ため、ヘシアンを計算している時間を使って構造最適化を進める方が有意義であることが良くあります。
そこで、最適化のトラジェクトリを使って短時間でヘシアンを近似するのが準ニュートン法です。
ヘシアンが厳密でないためNewton法ほど少ないステップ数で収束するというわけにはいきませんが、それでもかなり高速に収束します。
BFGS¶
BFGS法は準ニュートン法の中でも最も標準的に使われているヘシアン近似方法です。
他にも SR1法やDFP法等様々な方法が知られていますが、ASEではBFGS法が実装されています。
LBFGS¶
BFGS法の欠点として、最適化した変数の次元の大きさの行列を格納し行列ベクトル積を計算しないといけないというものがあります。
これは空間計算量が次元 \(N\) に対して \(O(N^2)\) 、時間計算量も \(O(N^2)\) となり、多次元の関数を最適化するのに不向きです。
そこで直近のいくつかのトラジェクトリの情報だけを用いてBFGS法とほとんど同じ結果になるように変形し、\(O(N)\) の計算量でBFGSをほとんど再現できるようにしたのがLBFGS法です。
PFPの計算量は \(O(N)\) なので、BFGSのような \(O(N^2)\) の手法を用いると原子数300程度からOptimizerでの計算時間が目立ってきますが、LBFGSだと原子数が増えてもOptimizerの計算時間は目立ってきません。
LineSearch¶
ポテンシャルが二次関数でない場合、ニュートン法は振動してしまうことがあります。また、準ニュートン法ではヘシアンが正確でないこともあり、適切なステップサイズを提示できないことがあります。
このような場合でも最適化を収束させるため、LineSearchを行うと最適化が安定します。
準ニュートン法で最適化の方向を決めた後、最適化のステップサイズに関してはその方向にエネルギーが十分小さくなるようなステップサイズをとることが望ましいです。このようなステップサイズを実際に何点か計算してみて決定するのがLineSearchです。
BFGS等だとニュートン法の方法に従って計算した結果エネルギーが上がってしまったというようなケースがありえますが、BFGSLineSearchだとそのようなケースでも良いステップサイズを探索します。
ベンチマーク¶
[1]:
import os
from time import perf_counter
from typing import Optional, Type
import matplotlib.pyplot as plt
import pandas as pd
from ase import Atoms
from ase.build import bulk, molecule
from ase.calculators.calculator import Calculator
from ase.io import Trajectory
from ase.optimize import BFGS, FIRE, LBFGS, BFGSLineSearch, LBFGSLineSearch, MDMin
from ase.optimize.optimize import Optimizer
from tqdm.auto import tqdm
import pfp_api_client
from pfcc_extras.structure.ase_rdkit_converter import atoms_to_smiles, smiles_to_atoms
from pfcc_extras.visualize.view import view_ngl
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
[2]:
calc = ASECalculator(Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v2.0.0"))
calc_mol = ASECalculator(Estimator(calc_mode=EstimatorCalcMode.MOLECULE, model_version="v2.0.0"))
print(f"pfp_api_client: {pfp_api_client.__version__}")
print(calc.estimator.model_version)
pfp_api_client: 1.18.0
v2.0.0
Force計算回数の比較¶
以下はCH3CHO分子の構造最適化を行う簡単な例です。
molecule("CH3CHO") で得られる構造は厳密な安定構造ではありませんが、かなり良い構造をしています。このような構造から開始した時には準ニュートン法を用いるとかなり少ない回数のPFP呼び出しで構造最適化を完了することが出来ます。
まずは構造最適化が完了するまでにForceを何回計算するかベンチマークをとってみましょう。
ここではPtのバルクを乱数で動かした構造を初期値として、最適構造になるまでのForce計算回数を測定します。
[3]:
def get_force_calls(opt: Optimizer) -> int:
"""Obtrain number of force calculations"""
if isinstance(opt, (BFGSLineSearch, LBFGSLineSearch)):
return opt.force_calls
else:
return opt.nsteps
[4]:
os.makedirs("output", exist_ok=True)
[5]:
atoms_0 = bulk("Pt") * (4, 4, 4)
atoms_0.rattle(stdev=0.1)
[6]:
view_ngl(atoms_0)
[7]:
force_calls = {}
for opt_class in tqdm([FIRE, BFGS, LBFGS, BFGSLineSearch, LBFGSLineSearch]):
atoms = atoms_0.copy()
atoms.calc = calc
name = opt_class.__name__
with opt_class(atoms, trajectory=f"output/{name}_ch3cho.traj", logfile=None) as opt:
opt.run(fmax=0.05)
force_calls[name] = [get_force_calls(opt)]
[8]:
df = pd.DataFrame.from_dict(force_calls)
df.plot.bar()
plt.ylabel("Force calls")
plt.show()
FIREが抜きんでてForce評価回数が多く、BFGSとLBFGSは同じくらいとなりました。
また、LineSearchのある手法もBFGSやLBFGSと同じくらいのForce評価回数となりました。
1ステップあたりの速度比較¶
次に同じそれぞれの最適化手法の1stepあたりの速度を比較してみます。
[9]:
from ase.build import bulk
atoms = bulk("Fe") * (10, 5, 3)
atoms.rattle(stdev=0.2)
atoms.calc = calc
view_ngl(atoms)
[10]:
def opt_benchmark(atoms: Atoms, calculator: Calculator, opt_class, logfile: Optional[str] = "-", steps: int = 10) -> float:
_atoms = atoms.copy()
_atoms.calc = calculator
with opt_class(_atoms, logfile=logfile) as opt:
start_time = perf_counter()
opt.run(steps=steps)
end_time = perf_counter()
duration = end_time - start_time
duration_per_force_calls = duration / get_force_calls(opt)
return duration_per_force_calls
[11]:
result_dict = {
"n_atoms": [],
"FIRE": [],
"BFGS": [],
"LBFGS": [],
"BFGSLineSearch": [],
"LBFGSLineSearch": [],
}
for i in range(1, 10):
atoms = bulk("Fe") * (10, 10, i)
atoms.rattle(stdev=0.2)
n_atoms = atoms.get_global_number_of_atoms()
result_dict["n_atoms"].append(n_atoms)
print(f"----- n_atoms {n_atoms} -----")
for opt_class in [FIRE, BFGS, LBFGS, BFGSLineSearch, LBFGSLineSearch]:
name = opt_class.__name__
print(f"Running {name}...")
duration = opt_benchmark(atoms, calc, opt_class=opt_class, logfile=None, steps=10)
print(f"Done in {duration:.2f} sec.")
result_dict[name].append(duration)
----- n_atoms 100 -----
Running FIRE...
Done in 0.07 sec.
Running BFGS...
Done in 0.08 sec.
Running LBFGS...
Done in 0.07 sec.
Running BFGSLineSearch...
Done in 0.07 sec.
Running LBFGSLineSearch...
Done in 0.07 sec.
----- n_atoms 200 -----
Running FIRE...
Done in 0.07 sec.
Running BFGS...
Done in 0.14 sec.
Running LBFGS...
Done in 0.07 sec.
Running BFGSLineSearch...
Done in 0.09 sec.
Running LBFGSLineSearch...
Done in 0.07 sec.
----- n_atoms 300 -----
Running FIRE...
Done in 0.08 sec.
Running BFGS...
Done in 0.25 sec.
Running LBFGS...
Done in 0.13 sec.
Running BFGSLineSearch...
Done in 0.14 sec.
Running LBFGSLineSearch...
Done in 0.07 sec.
----- n_atoms 400 -----
Running FIRE...
Done in 0.09 sec.
Running BFGS...
Done in 0.45 sec.
Running LBFGS...
Done in 0.09 sec.
Running BFGSLineSearch...
Done in 0.24 sec.
Running LBFGSLineSearch...
Done in 0.08 sec.
----- n_atoms 500 -----
Running FIRE...
Done in 0.10 sec.
Running BFGS...
Done in 0.77 sec.
Running LBFGS...
Done in 0.10 sec.
Running BFGSLineSearch...
Done in 0.33 sec.
Running LBFGSLineSearch...
Done in 0.09 sec.
----- n_atoms 600 -----
Running FIRE...
Done in 0.11 sec.
Running BFGS...
Done in 1.37 sec.
Running LBFGS...
Done in 0.11 sec.
Running BFGSLineSearch...
Done in 0.46 sec.
Running LBFGSLineSearch...
Done in 0.11 sec.
----- n_atoms 700 -----
Running FIRE...
Done in 0.13 sec.
Running BFGS...
Done in 1.98 sec.
Running LBFGS...
Done in 0.11 sec.
Running BFGSLineSearch...
Done in 0.68 sec.
Running LBFGSLineSearch...
Done in 0.11 sec.
----- n_atoms 800 -----
Running FIRE...
Done in 0.15 sec.
Running BFGS...
Done in 2.96 sec.
Running LBFGS...
Done in 0.14 sec.
Running BFGSLineSearch...
Done in 0.88 sec.
Running LBFGSLineSearch...
Done in 0.12 sec.
----- n_atoms 900 -----
Running FIRE...
Done in 0.15 sec.
Running BFGS...
Done in 4.07 sec.
Running LBFGS...
Done in 0.14 sec.
Running BFGSLineSearch...
Done in 1.07 sec.
Running LBFGSLineSearch...
Done in 0.14 sec.
[12]:
import matplotlib.pyplot as plt
import pandas as pd
fig, axes = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))
df = pd.DataFrame(result_dict)
df = df.set_index("n_atoms")
df.plot(ax=axes[0])
axes[0].set_title("Speed for various opt algorithms")
df[["LBFGS", "FIRE", "LBFGSLineSearch"]].plot(ax=axes[1])
axes[1].set_title("Speed for various opt algorithms ($O(N)$ methods only)")
for ax in axes:
ax.set_ylabel("sec / step")
fig.savefig("output/opt_benchmark.png")
plt.show(fig)
BFGSやBFGSLineSearchでは原子数が多くなるにつれて1ステップあたりの計算時間が長くなってしまっています。
ASEに実装されているBFGS法ではHessianの逆行列を求める際に固有値分解を毎ステップごとに行います。
Hessianは縦横それぞれ原子数をNとして3Nの大きさをもっており、固有値分解の計算量は \(O(N^3)\) であるため、ASEのBFGS法は原子数の3乗のオーダーで急速に遅くなっていきます。
また、BFGSLineSearchではHessianの固有値分解を行うことなく、Hessianの逆行列を直接近似します。このため、行列ベクトル積が律速になり \(O(N^2)\) とBFGSよりは速くなります。
しかし、原子数が多い時には \(O(N)\) で動作するLBFGSやFIREよりはずっと遅くなります。
最適化アルゴリズムの使い分け¶
これまでに、FIREは他の手法よりもForce呼び出し回数が多く、BFGSやBFGSLineSearchは原子数が多くなるとForceの評価よりもOptimizerの内部ルーチンの方が律速になることを見てきました。
ではどういった場合でもLBFGSとLBFGSLineSearchを用いると良いのでしょうか？
これからそれらのアルゴリズムでもうまくいかない例もあることを見ていきましょう。
LBFGSやBFGSがうまくいかない例¶
次の例を見てみましょう。 この例ではとても単純なポテンシャルである Lennard-Jones (LJ) potentialを用いて、その安定点から離れたところを初期値として構造最適化を行っています。
LJ potential は次式で表される形をしており、以下のような形をしています。
\[\phi (r) = 4 \epsilon \bigg\{ \left( \frac{\sigma}{r} \right) ^{12} - \left( \frac{\sigma}{r} \right) ^{6} \bigg\}\]
[13]:
from ase.calculators.lj import LennardJones
import numpy as np
calc_lj = LennardJones()
dists = np.linspace(1.0, 4.0, 500)
E_pot_list = []
for d in dists:
atoms = Atoms("Ar2", [[0, 0, 0], [d, 0, 0]])
atoms.calc = calc_lj
E_pot = atoms.get_potential_energy()
E_pot_list.append(E_pot)
plt.plot(dists, E_pot_list)
plt.title("Lennard Jones Potential")
plt.ylabel("eV")
plt.xlabel("Distance [A]")
plt.show()
この単純なPotential の形でDistanceを安定点である1.12A から離れたところからスタートし構造最適化をしてみましょう。
[14]:
def lennard_jones_trajectory(opt_class: Type[Optimizer], d0: float):
name = opt_class.__name__
trajectory_path = f"output/Ar_{name}.traj"
atoms = Atoms("Ar2", [[0, 0, 0], [d0, 0, 0]])
atoms.calc = calc_lj
opt = opt_class(atoms, trajectory=trajectory_path)
opt.run()
distance_list = []
energy_list = []
for atoms in Trajectory(trajectory_path):
energy_list.append(atoms.get_potential_energy())
distance_list.append(atoms.get_distance(0, 1))
print("Distance in opt trajectory: ", distance_list)
plt.plot(dists, E_pot_list)
plt.plot(distance_list, energy_list, marker="x")
plt.title(f"Lennard Jones Potential: {name} opt")
plt.ylim(-1.0, 0)
plt.ylabel("eV")
plt.xlabel("Distance [A]")
plt.show()
まずはFIREを行った場合です。安定して最適値1.12A にたどり着くことが出来ています。
[15]:
lennard_jones_trajectory(FIRE, 1.5)
Step Time Energy fmax
FIRE: 0 14:15:06 -0.314857 1.158029
FIRE: 1 14:15:06 -0.342894 1.264380
FIRE: 2 14:15:06 -0.410024 1.512419
FIRE: 3 14:15:06 -0.546744 1.968308
FIRE: 4 14:15:06 -0.812180 2.383891
FIRE: 5 14:15:06 -0.862167 5.585844
FIRE: 6 14:15:06 -0.966377 2.149190
FIRE: 7 14:15:06 -0.991817 0.522302
FIRE: 8 14:15:06 -0.992148 0.491275
FIRE: 9 14:15:06 -0.992732 0.429940
FIRE: 10 14:15:06 -0.993427 0.339973
FIRE: 11 14:15:06 -0.994057 0.224455
FIRE: 12 14:15:06 -0.994451 0.088426
FIRE: 13 14:15:06 -0.994489 0.060610
FIRE: 14 14:15:06 -0.994490 0.059506
FIRE: 15 14:15:07 -0.994492 0.057320
FIRE: 16 14:15:07 -0.994495 0.054094
FIRE: 17 14:15:07 -0.994499 0.049889
Distance in opt trajectory: [1.5, 1.4768394233790767, 1.42839123774426, 1.349694678243147, 1.231631957294667, 1.0658914242047832, 1.0938206432760949, 1.132495812583008, 1.1318429355798734, 1.1305759647542164, 1.1287715691339621, 1.1265422074054834, 1.1240322773166505, 1.1214118150599879, 1.1214307556785812, 1.1214682920288386, 1.1215237409768002, 1.1215960942672192]
[16]:
lennard_jones_trajectory(FIRE, 1.28)
Step Time Energy fmax
FIRE: 0 14:15:07 -0.697220 2.324552
FIRE: 1 14:15:07 -0.807702 2.387400
FIRE: 2 14:15:07 -0.987241 0.821987
FIRE: 3 14:15:07 -0.520102 13.570402
FIRE: 4 14:15:07 -0.971689 1.903593
FIRE: 5 14:15:07 -0.939109 1.840004
FIRE: 6 14:15:07 -0.943288 1.793714
FIRE: 7 14:15:07 -0.951214 1.694594
FIRE: 8 14:15:07 -0.961964 1.529352
FIRE: 9 14:15:07 -0.974030 1.277952
FIRE: 10 14:15:07 -0.985244 0.914936
FIRE: 11 14:15:07 -0.992867 0.414203
FIRE: 12 14:15:07 -0.994043 0.239611
FIRE: 13 14:15:07 -0.994061 0.235001
FIRE: 14 14:15:07 -0.994095 0.225890
FIRE: 15 14:15:07 -0.994143 0.212487
FIRE: 16 14:15:07 -0.994201 0.195100
FIRE: 17 14:15:07 -0.994265 0.174121
FIRE: 18 14:15:07 -0.994330 0.150013
FIRE: 19 14:15:07 -0.994391 0.123296
FIRE: 20 14:15:07 -0.994449 0.091418
FIRE: 21 14:15:07 -0.994495 0.054286
FIRE: 22 14:15:07 -0.994519 0.012221
Distance in opt trajectory: [1.28, 1.2335089629222111, 1.1392699328939198, 1.0285911631178963, 1.0964431714345662, 1.1738131464936081, 1.1715131413850235, 1.1669709941864936, 1.160310604717144, 1.1517385248429906, 1.1415690055780352, 1.130255815977825, 1.118424872251103, 1.1184997505580012, 1.1186480666965142, 1.1188669733107373, 1.1191522821112356, 1.1194985597723826, 1.1198992503739216, 1.1203468201911215, 1.1208857684008715, 1.121520438052656, 1.1222486280893154]
BFGS法を1.5Aから開始した場合、特定のstepではあまりにも大きくジャンプしてしまい、安定点を大きく飛び越えるような不安定な挙動になっていることがわかります。 (1.18A –> 0.78Aにジャンプしているようなところです。) このように、最適値1.12A にたどり着くまでに、この単純な例でも大きく振動しまいます。 これは2次の近似を行う手法では、関数曲面が2次関数と異なる場合に、その極小値の見積もりが大きくはずれてしまうために起こります。
[17]:
lennard_jones_trajectory(BFGS, 1.5)
Step Time Energy fmax
BFGS: 0 14:15:07 -0.314857 1.158029
BFGS: 1 14:15:07 -0.355681 1.312428
BFGS: 2 14:15:07 -0.916039 2.041011
BFGS: 3 14:15:07 55.304434 974.487533
BFGS: 4 14:15:07 -0.917740 2.028761
BFGS: 5 14:15:07 -0.919418 2.016337
BFGS: 6 14:15:07 -0.747080 8.517609
BFGS: 7 14:15:07 -0.965068 1.472894
BFGS: 8 14:15:07 -0.984660 0.939655
BFGS: 9 14:15:07 -0.992324 0.528593
BFGS: 10 14:15:07 -0.994444 0.092418
BFGS: 11 14:15:07 -0.994520 0.007274
Distance in opt trajectory: [1.5, 1.4669134619701096, 1.1856710148074785, 0.7856710148074794, 1.1848349876531552, 1.1840057048694352, 1.0494144415981015, 1.1582431525834522, 1.1421985974452062, 1.1139254943723917, 1.1241042745909042, 1.1225894828343739]
同じことがLBFGS法でも起こります。
このような振動が多原子系で起こった場合、本当にいつまでたっても最適化が終了しないような場合もあります。
[18]:
lennard_jones_trajectory(LBFGS, 1.28)
Step Time Energy fmax
LBFGS: 0 14:15:07 -0.697220 2.324552
LBFGS: 1 14:15:07 -0.854683 2.314904
LBFGS: 2 14:15:07 33.771592 599.751912
LBFGS: 3 14:15:07 -0.858237 2.305680
LBFGS: 4 14:15:07 -0.861748 2.295944
LBFGS: 5 14:15:07 17.305333 318.491824
LBFGS: 6 14:15:07 -0.867638 2.278148
LBFGS: 7 14:15:07 -0.873395 2.258876
LBFGS: 8 14:15:07 5.604073 121.154807
LBFGS: 9 14:15:07 -0.885566 2.211348
LBFGS: 10 14:15:07 -0.897001 2.157091
LBFGS: 11 14:15:07 0.369118 30.787165
LBFGS: 12 14:15:07 -0.925197 1.970726
LBFGS: 13 14:15:07 -0.947015 1.749104
LBFGS: 14 14:15:07 -0.915603 4.009032
LBFGS: 15 14:15:07 -0.985429 0.906882
LBFGS: 16 14:15:07 -0.993163 0.377052
LBFGS: 17 14:15:07 -0.994465 0.080647
LBFGS: 18 14:15:07 -0.994520 0.005364
Distance in opt trajectory: [1.28, 1.213584232746016, 0.813584232746016, 1.2120462614888505, 1.2105202846657441, 0.8506725019776178, 1.2079447803538972, 1.2054073889508412, 0.9079986088524912, 1.1999638317862706, 1.1947303311793407, 0.9866659932915066, 1.1811069055371053, 1.169409273526655, 1.0770879993440852, 1.1413655567592533, 1.1295077111909069, 1.1210690913427286, 1.1225559889249785]
また、他にも初期値が安定点から離れたところからスタートした場合に正しく安定点にたどり着けていないケースも有ります。 これは、Hessianの初期化がうまく行かないなどの理由で起きていると考えられます。
[19]:
lennard_jones_trajectory(LBFGS, 1.50)
Step Time Energy fmax
LBFGS: 0 14:15:08 -0.314857 1.158029
LBFGS: 1 14:15:08 -0.355681 1.312428
LBFGS: 2 14:15:08 -0.129756 0.447300
LBFGS: 3 14:15:08 -0.079409 0.263016
LBFGS: 4 14:15:08 -0.040472 0.129677
LBFGS: 5 14:15:08 -0.021155 0.068924
LBFGS: 6 14:15:08 -0.009646 0.035707
Distance in opt trajectory: [1.5, 1.4669134619701096, 1.748155909132741, 1.8935676545167726, 2.101103976550205, 2.302941039822042, 2.531922335240822]
このようにBFGSやLBFGSがうまくいかない例であっても、LineSearchを行うことで正しく収束させることが出来ます。
LBFGSLineSearchはほとんどの場合において最良のOptimizerです。
[20]:
lennard_jones_trajectory(LBFGSLineSearch, 1.5)
Step Time Energy fmax
LBFGSLineSearch: 0 14:15:08 -0.314857 1.158029
LBFGSLineSearch: 1 14:15:08 -0.994517 0.021446
Distance in opt trajectory: [1.5, 1.1220880716401869]
LBFGSLineSearchでは1ステップで終了しているように見えますが、内部でLineSearchを行っている部分がステップにカウントされていないため、実際の計算数はもう少し多いです。 BFGSLineSearchでは Step[FC]という表示がありますが、このFCが実際のフォース評価回数です。
[21]:
lennard_jones_trajectory(BFGSLineSearch, 1.5)
Step[ FC] Time Energy fmax
BFGSLineSearch: 0[ 0] 14:15:08 -0.314857 1.1580
BFGSLineSearch: 1[ 5] 14:15:08 -0.994057 0.2244
BFGSLineSearch: 2[ 7] 14:15:08 -0.994520 0.0054
Distance in opt trajectory: [1.5, 1.1265420003736115, 1.1225561293552997]
ここではBFGSやLBFGSがうまくいかない例があることを見てきました。 しかし、このような例であってもLBFGSLineSearchはうまくいっています。
ここまでの結論としては多くの場合においてLBFGSLineSearchが最良のOptimizerとなります。
以下にLBFGSLineSearch以外のOptimizerの欠点を述べます。
FIRE: 収束までのステップ数が長い
BFGS: 原子数が多い時に計算時間がかかる。
BFGSLineSearch: 原子数が多い時に計算時間がかかる。
LBFGS: 振動したり、エネルギーが上がってしまう場合がある。
しかし、常にLBFGSLineSearchを使うとうまくいくかというと、現状残念ながら必ずしもそうではありません。
LBFGSLineSearchが良くない例¶
LBFGSLineSearchは非常に完成度の高い手法であり、多くの場合はLBFGSLineSearchを使えばうまくいきます。
しかし、稀にLBFGSLineSearchがうまくいかなくなってしまう例も存在します。
例えば、fmaxを小さくして精密に構造最適化したい時に良くこのような例にあたります。
例えばトルエンのメチル基の回転を最適化したい場合、fmax=0.001eV/A程度にしないと回転に関して極小値に到達しません。
この例ではLineSearchがなかなか収束しなくなり1ステップに長い時間がかかるようになったり、エネルギーが下がらなくなったり、場合によってはLineSearchが破綻して
RuntimeError: LineSearch failed! となってエラーになってしまう場合もあります。[22]:
atoms_0 = smiles_to_atoms("Cc1ccccc1")
tmp = atoms_0[7:10]
tmp.rotate([1.0, 0.0, 0.0], 15.0)
atoms_0.positions[7:10] = tmp.positions
atoms_0.calc = calc_mol
[23]:
view_ngl(atoms_0, representations=["ball+stick"], h=600, w=600)
[24]:
# Please set '-', if you want to see detailed logs.
logfile = None
steps = {}
images = []
for optimizer_class in (FIRE, LBFGS, LBFGSLineSearch):
name = optimizer_class.__name__
atoms = atoms_0.copy()
atoms.calc = calc_mol
with optimizer_class(atoms, logfile=logfile) as opt:
try:
print(f"{name} optimization starts.")
opt.run(fmax=0.001, steps=200)
print(f"Optimization finished without error. steps = {opt.nsteps}")
finally:
steps[name] = [get_force_calls(opt)]
images.append(atoms.copy())
FIRE optimization starts.
Optimization finished without error. steps = 200
LBFGS optimization starts.
Optimization finished without error. steps = 80
LBFGSLineSearch optimization starts.
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/tmp/ipykernel_1863/983813568.py in <cell line: 6>()
11 try:
12 print(f"{name} optimization starts.")
---> 13 opt.run(fmax=0.001, steps=200)
14 print(f"Optimization finished without error. steps = {opt.nsteps}")
15 finally:
~/.py39/lib/python3.9/site-packages/ase/optimize/optimize.py in run(self, fmax, steps)
428 """
429 self.fmax = fmax
--> 430 return Dynamics.run(self, steps=steps)
431
432 def converged(self, forces=None):
~/.py39/lib/python3.9/site-packages/ase/optimize/optimize.py in run(self, steps)
273 """
274
--> 275 for converged in Dynamics.irun(self, steps=steps):
276 pass
277 return converged
~/.py39/lib/python3.9/site-packages/ase/optimize/optimize.py in irun(self, steps)
244 while not is_converged and self.nsteps < self.max_steps:
245 # compute the next step
--> 246 self.step()
247 self.nsteps += 1
248
~/.py39/lib/python3.9/site-packages/ase/optimize/lbfgs.py in step(self, forces)
154 if self.use_line_search is True:
155 e = self.func(pos)
--> 156 self.line_search(pos, g, e)
157 dr = (self.alpha_k * self.p).reshape(len(self.optimizable), -1)
158 else:
~/.py39/lib/python3.9/site-packages/ase/optimize/lbfgs.py in line_search(self, r, g, e)
247 c2=.46, stpmax=50.)
248 if self.alpha_k is None:
--> 249 raise RuntimeError('LineSearch failed!')
250
251
RuntimeError: LineSearch failed!
[25]:
steps["index"] = "nsteps"
df = pd.DataFrame(steps)
df = df.set_index("index")
df
[25]:
| FIRE | LBFGS | LBFGSLineSearch | |
|---|---|---|---|
| index | |||
| nsteps | 200 | 80 | 439 |
この例ではLBFGSは収束しますが、FIREは収束しませんでした。
また、LBFGSLineSearchは後半に1ステップあたりにすごく長い時間がかかった上に収束せず、場合によっては
RuntimeError: LineSearch failed! というエラーが出ます。経験的に多くの場合はfmax=0.002eV/A 付近からLBFGSLineSearchが収束しにくくなりますが、
もっとポテンシャルの形状が悪い場合はfmax=0.05に到達しないケースもあります。
ちなみに、このようにLBFGSLineSearchが進まなくなってしまった場合はLBFGSLineSearchを実行しなおすことで、うまくいくようになることがあります。（うまくいかない時もあります）
[26]:
with LBFGSLineSearch(atoms) as opt:
opt.run(0.001)
Step Time Energy fmax
LBFGSLineSearch: 0 14:16:21 -72.184885 0.046849
LBFGSLineSearch: 1 14:16:21 -72.184900 0.009447
LBFGSLineSearch: 2 14:16:21 -72.184905 0.006908
LBFGSLineSearch: 3 14:16:21 -72.184910 0.006823
LBFGSLineSearch: 4 14:16:21 -72.184913 0.003909
LBFGSLineSearch: 5 14:16:22 -72.184913 0.003907
LBFGSLineSearch: 6 14:16:22 -72.184916 0.005214
LBFGSLineSearch: 7 14:16:22 -72.184918 0.003444
LBFGSLineSearch: 8 14:16:23 -72.184920 0.003443
LBFGSLineSearch: 9 14:16:24 -72.184923 0.003432
LBFGSLineSearch: 10 14:16:24 -72.184926 0.003433
LBFGSLineSearch: 11 14:16:24 -72.184925 0.003432
LBFGSLineSearch: 12 14:16:25 -72.184926 0.003422
LBFGSLineSearch: 13 14:16:25 -72.184925 0.003421
LBFGSLineSearch: 14 14:16:26 -72.184925 0.003422
LBFGSLineSearch: 15 14:16:26 -72.184925 0.003422
LBFGSLineSearch: 16 14:16:26 -72.184924 0.003423
LBFGSLineSearch: 17 14:16:27 -72.184924 0.003421
LBFGSLineSearch: 18 14:16:27 -72.184924 0.003421
LBFGSLineSearch: 19 14:16:28 -72.184925 0.003421
LBFGSLineSearch: 20 14:16:28 -72.184924 0.003421
LBFGSLineSearch: 21 14:16:28 -72.184922 0.003420
LBFGSLineSearch: 22 14:16:29 -72.184926 0.003420
LBFGSLineSearch: 23 14:16:29 -72.184926 0.003421
LBFGSLineSearch: 24 14:16:29 -72.184925 0.003421
LBFGSLineSearch: 25 14:16:29 -72.184924 0.003421
LBFGSLineSearch: 26 14:16:30 -72.184927 0.003421
LBFGSLineSearch: 27 14:16:30 -72.184925 0.003421
LBFGSLineSearch: 28 14:16:31 -72.184924 0.003420
LBFGSLineSearch: 29 14:16:31 -72.184926 0.003421
LBFGSLineSearch: 30 14:16:31 -72.184925 0.003420
LBFGSLineSearch: 31 14:16:32 -72.184923 0.003421
LBFGSLineSearch: 32 14:16:32 -72.184926 0.003421
LBFGSLineSearch: 33 14:16:32 -72.184924 0.003420
LBFGSLineSearch: 34 14:16:33 -72.184922 0.003421
LBFGSLineSearch: 35 14:16:33 -72.184925 0.003421
LBFGSLineSearch: 36 14:16:33 -72.184921 0.003421
LBFGSLineSearch: 37 14:16:34 -72.184924 0.003420
LBFGSLineSearch: 38 14:16:34 -72.184924 0.003421
LBFGSLineSearch: 39 14:16:34 -72.184923 0.003421
LBFGSLineSearch: 40 14:16:34 -72.184921 0.003421
LBFGSLineSearch: 41 14:16:34 -72.184926 0.003421
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/tmp/ipykernel_1863/1423403903.py in <cell line: 1>()
1 with LBFGSLineSearch(atoms) as opt:
----> 2 opt.run(0.001)
~/.py39/lib/python3.9/site-packages/ase/optimize/optimize.py in run(self, fmax, steps)
428 """
429 self.fmax = fmax
--> 430 return Dynamics.run(self, steps=steps)
431
432 def converged(self, forces=None):
~/.py39/lib/python3.9/site-packages/ase/optimize/optimize.py in run(self, steps)
273 """
274
--> 275 for converged in Dynamics.irun(self, steps=steps):
276 pass
277 return converged
~/.py39/lib/python3.9/site-packages/ase/optimize/optimize.py in irun(self, steps)
244 while not is_converged and self.nsteps < self.max_steps:
245 # compute the next step
--> 246 self.step()
247 self.nsteps += 1
248
~/.py39/lib/python3.9/site-packages/ase/optimize/lbfgs.py in step(self, forces)
154 if self.use_line_search is True:
155 e = self.func(pos)
--> 156 self.line_search(pos, g, e)
157 dr = (self.alpha_k * self.p).reshape(len(self.optimizable), -1)
158 else:
~/.py39/lib/python3.9/site-packages/ase/optimize/lbfgs.py in line_search(self, r, g, e)
247 c2=.46, stpmax=50.)
248 if self.alpha_k is None:
--> 249 raise RuntimeError('LineSearch failed!')
250
251
RuntimeError: LineSearch failed!
[27]:
view_ngl(images, representations=["ball+stick"])
上記知見を踏まえた、最適化アルゴリズムの実装¶
LBFGSLineSearchは多くの場合最良のOptimizerですが、fmaxが小さいと非常に遅くなってしまうだけでなく、場合によってはエラーが出てしまうことがあります。
FIREはLBFGSLineSearchと比べて遅いことが多いですが、多くの場合頑健に動作します。ただし、fmaxが小さいとなかなか収束しないことがあります。
LBFGSはfmaxが小さくても収束しますが、Lennard-Jonesの例で見たように振動したりエネルギーが上がってしまうケースがあります。
ヒューリスティックな表現をすると、2次関数近似が良い近似となる安定点に近いところでは、収束性の良い2次のニュートン法を用いたいです。 ただし、LBFGSLineSearchが数値的に安定しないケースがあるためLBFGSを使いたいです。
２次関数近似がよい近似とならない不安定な構造では、FIREやLBFGSLineSearchを使って頑健にエネルギーを落としたいです。
しかし、LBFGSLineSearchが数値的に安定しないケースがどれくらのfmaxで起こるか予想しきれないため、 ハイスループット計算など、未知の物質に対して統一的に手法を適用する場合ではFIREのほうが安定に動作します。
ここまでの知見をまとめると、ヒューリスティックにうまくいきそうな最適化アルゴリズムとして、
構造最適化の初期など、不安定な場所では、FIREを用いて勾配を下っていく。
ある程度安定点に達したら(Forceなどが小さな値になったら)LBFGS法を用いて高速に収束を試みる。
という方法がうまくいきそうです。
matlantis-features では、このようなアルゴリズムを FIRELBFGS Optimizer として実装をしています。
[28]:
from matlantis_features.ase_ext.optimize import FIRELBFGS
まずはLBFGSLineSearchが破綻した例であるトルエンのメチル基の回転を見てみましょう。
このような例でも破綻することなく、LBFGSよりは遅いがFIREよりはずっと速く、頑健に収束しています。
[29]:
atoms_0 = smiles_to_atoms("Cc1ccccc1")
tmp = atoms_0[7:10]
tmp.rotate([1.0, 0.0, 0.0], 15.0)
atoms_0.positions[7:10] = tmp.positions
atoms_0.calc = calc_mol
[30]:
view_ngl(atoms_0, representations=["ball+stick"], h=600, w=600)
[31]:
steps = {}
images = []
atoms = atoms_0.copy()
atoms.calc = calc_mol
with FIRELBFGS(atoms, logfile="-") as opt:
try:
opt.run(fmax=0.001, steps=200)
finally:
steps[name] = [get_force_calls(opt)]
images.append(atoms.copy())
Step Time Energy fmax
FIRELBFGS: 0 14:16:52 -72.042224 1.929766
FIRELBFGS: 1 14:16:52 -72.115609 0.631807
FIRELBFGS: 2 14:16:52 -72.123116 1.466699
FIRELBFGS: 3 14:16:52 -72.136854 1.138275
FIRELBFGS: 4 14:16:52 -72.153085 0.575332
FIRELBFGS: 5 14:16:52 -72.160511 0.396632
FIRELBFGS: 6 14:16:53 -72.159784 0.626498
FIRELBFGS: 7 14:16:53 -72.160875 0.593301
FIRELBFGS: 8 14:16:53 -72.162798 0.528663
FIRELBFGS: 9 14:16:53 -72.165194 0.436117
FIRELBFGS: 10 14:16:53 -72.167603 0.321110
FIRELBFGS: 11 14:16:53 -72.169658 0.225238
FIRELBFGS: 12 14:16:53 -72.171126 0.206753
FIRELBFGS: 13 14:16:53 -72.172131 0.228274
FIRELBFGS: 14 14:16:53 -72.172988 0.253087
FIRELBFGS: 15 14:16:53 -72.174026 0.302716
FIRELBFGS: 16 14:16:53 -72.175554 0.345937
FIRELBFGS: 17 14:16:53 -72.177558 0.327525
FIRELBFGS: 18 14:16:53 -72.179640 0.236171
FIRELBFGS: 19 14:16:53 -72.181110 0.167093
FIRELBFGS: 20 14:16:53 -72.181669 0.174786
FIRELBFGS: 21 14:16:53 -72.181748 0.162253
FIRELBFGS: 22 14:16:53 -72.181878 0.138127
FIRELBFGS: 23 14:16:53 -72.182052 0.104237
FIRELBFGS: 24 14:16:53 -72.182208 0.069203
FIRELBFGS: 25 14:16:54 -72.182334 0.066081
FIRELBFGS: 26 14:16:54 -72.182405 0.062471
FIRELBFGS: 27 14:16:54 -72.182446 0.068936
FIRELBFGS: 28 14:16:54 -72.182486 0.093038
FIRELBFGS: 29 14:16:54 -72.182539 0.102686
FIRELBFGS: 30 14:16:54 -72.182651 0.095438
FIRELBFGS: 31 14:16:54 -72.182820 0.073343
FIRELBFGS: 32 14:16:54 -72.183015 0.049743
FIRELBFGS: 33 14:16:54 -72.183170 0.040306
FIRELBFGS: 34 14:16:54 -72.183288 0.071183
FIRELBFGS: 35 14:16:54 -72.183422 0.072146
FIRELBFGS: 36 14:16:54 -72.183582 0.048024
FIRELBFGS: 37 14:16:54 -72.183729 0.039146
FIRELBFGS: 38 14:16:54 -72.183804 0.054861
FIRELBFGS: 39 14:16:54 -72.183881 0.046277
FIRELBFGS: 40 14:16:54 -72.183929 0.049557
FIRELBFGS: 41 14:16:54 -72.183948 0.040189
FIRELBFGS: 42 14:16:54 -72.183970 0.026913
FIRELBFGS: 43 14:16:55 -72.183991 0.025744
FIRELBFGS: 44 14:16:55 -72.184008 0.024797
FIRELBFGS: 45 14:16:55 -72.184019 0.031188
FIRELBFGS: 46 14:16:55 -72.184041 0.030759
FIRELBFGS: 47 14:16:55 -72.184068 0.022572
FIRELBFGS: 48 14:16:55 -72.184099 0.020500
FIRELBFGS: 49 14:16:55 -72.184124 0.017614
FIRELBFGS: 50 14:16:55 -72.184154 0.015145
FIRELBFGS: 51 14:16:55 -72.184190 0.019245
FIRELBFGS: 52 14:16:55 -72.184386 0.029275
FIRELBFGS: 53 14:16:55 -72.184437 0.024904
FIRELBFGS: 54 14:16:55 -72.184492 0.027085
FIRELBFGS: 55 14:16:55 -72.184535 0.022539
FIRELBFGS: 56 14:16:55 -72.184618 0.035791
FIRELBFGS: 57 14:16:55 -72.184682 0.025345
FIRELBFGS: 58 14:16:55 -72.184725 0.016344
FIRELBFGS: 59 14:16:55 -72.184730 0.011044
FIRELBFGS: 60 14:16:56 -72.184740 0.007274
FIRELBFGS: 61 14:16:56 -72.184745 0.006957
FIRELBFGS: 62 14:16:56 -72.184765 0.007375
FIRELBFGS: 63 14:16:56 -72.184769 0.010790
FIRELBFGS: 64 14:16:56 -72.184784 0.009786
FIRELBFGS: 65 14:16:56 -72.184792 0.007960
FIRELBFGS: 66 14:16:56 -72.184803 0.008439
FIRELBFGS: 67 14:16:56 -72.184809 0.010082
FIRELBFGS: 68 14:16:56 -72.184820 0.009707
FIRELBFGS: 69 14:16:56 -72.184823 0.005612
FIRELBFGS: 70 14:16:56 -72.184826 0.003727
FIRELBFGS: 71 14:16:56 -72.184825 0.003739
FIRELBFGS: 72 14:16:56 -72.184837 0.004489
FIRELBFGS: 73 14:16:56 -72.184836 0.004152
FIRELBFGS: 74 14:16:56 -72.184841 0.003730
FIRELBFGS: 75 14:16:56 -72.184833 0.003313
FIRELBFGS: 76 14:16:56 -72.184836 0.004958
FIRELBFGS: 77 14:16:56 -72.184835 0.005354
FIRELBFGS: 78 14:16:57 -72.184839 0.004680
FIRELBFGS: 79 14:16:57 -72.184843 0.004109
FIRELBFGS: 80 14:16:57 -72.184844 0.005151
FIRELBFGS: 81 14:16:57 -72.184850 0.005093
FIRELBFGS: 82 14:16:57 -72.184845 0.004273
FIRELBFGS: 83 14:16:57 -72.184853 0.004361
FIRELBFGS: 84 14:16:57 -72.184863 0.005749
FIRELBFGS: 85 14:16:57 -72.184860 0.009415
FIRELBFGS: 86 14:16:57 -72.184870 0.012794
FIRELBFGS: 87 14:16:57 -72.184886 0.012712
FIRELBFGS: 88 14:16:57 -72.184905 0.013424
FIRELBFGS: 89 14:16:57 -72.184930 0.013748
FIRELBFGS: 90 14:16:57 -72.184959 0.020716
FIRELBFGS: 91 14:16:57 -72.185021 0.032262
FIRELBFGS: 92 14:16:57 -72.185104 0.038434
FIRELBFGS: 93 14:16:57 -72.185215 0.038006
FIRELBFGS: 94 14:16:57 -72.185375 0.040919
FIRELBFGS: 95 14:16:57 -72.185574 0.043692
FIRELBFGS: 96 14:16:58 -72.184194 0.212188
FIRELBFGS: 97 14:16:58 -72.185368 0.128039
FIRELBFGS: 98 14:16:58 -72.185919 0.149830
FIRELBFGS: 99 14:16:58 -72.185997 0.128264
FIRELBFGS: 100 14:16:58 -72.186109 0.088823
FIRELBFGS: 101 14:16:58 -72.186193 0.038958
FIRELBFGS: 102 14:16:58 -72.186240 0.034949
FIRELBFGS: 103 14:16:58 -72.186253 0.050910
FIRELBFGS: 104 14:16:58 -72.186261 0.048247
FIRELBFGS: 105 14:16:58 -72.186268 0.043067
FIRELBFGS: 106 14:16:58 -72.186281 0.035627
FIRELBFGS: 107 14:16:58 -72.186298 0.026843
FIRELBFGS: 108 14:16:58 -72.186305 0.023197
FIRELBFGS: 109 14:16:58 -72.186314 0.019699
FIRELBFGS: 110 14:16:58 -72.186314 0.016956
FIRELBFGS: 111 14:16:58 -72.186331 0.017175
FIRELBFGS: 112 14:16:58 -72.186333 0.025444
FIRELBFGS: 113 14:16:58 -72.186338 0.030209
FIRELBFGS: 114 14:16:58 -72.186345 0.030242
FIRELBFGS: 115 14:16:58 -72.186347 0.024685
FIRELBFGS: 116 14:16:59 -72.186367 0.008148
FIRELBFGS: 117 14:16:59 -72.186368 0.007604
FIRELBFGS: 118 14:16:59 -72.186384 0.016938
FIRELBFGS: 119 14:16:59 -72.186391 0.014594
FIRELBFGS: 120 14:16:59 -72.186399 0.006488
FIRELBFGS: 121 14:16:59 -72.186412 0.005894
FIRELBFGS: 122 14:16:59 -72.186406 0.006581
FIRELBFGS: 123 14:16:59 -72.186417 0.007906
FIRELBFGS: 124 14:16:59 -72.186408 0.005107
FIRELBFGS: 125 14:16:59 -72.186420 0.003437
FIRELBFGS: 126 14:16:59 -72.186420 0.004048
FIRELBFGS: 127 14:16:59 -72.186423 0.005805
FIRELBFGS: 128 14:16:59 -72.186424 0.004949
FIRELBFGS: 129 14:16:59 -72.186428 0.003298
FIRELBFGS: 130 14:16:59 -72.186432 0.002626
FIRELBFGS: 131 14:16:59 -72.186425 0.003264
FIRELBFGS: 132 14:16:59 -72.186433 0.003322
FIRELBFGS: 133 14:16:59 -72.186431 0.003426
FIRELBFGS: 134 14:16:59 -72.186431 0.003106
FIRELBFGS: 135 14:17:00 -72.186437 0.003114
FIRELBFGS: 136 14:17:00 -72.186432 0.004069
FIRELBFGS: 137 14:17:00 -72.186439 0.005096
FIRELBFGS: 138 14:17:00 -72.186444 0.005440
FIRELBFGS: 139 14:17:00 -72.186440 0.005443
FIRELBFGS: 140 14:17:00 -72.186443 0.007676
FIRELBFGS: 141 14:17:00 -72.186451 0.008253
FIRELBFGS: 142 14:17:00 -72.186455 0.005408
FIRELBFGS: 143 14:17:00 -72.186465 0.005941
FIRELBFGS: 144 14:17:00 -72.186462 0.006734
FIRELBFGS: 145 14:17:00 -72.186461 0.011040
FIRELBFGS: 146 14:17:00 -72.186471 0.010580
FIRELBFGS: 147 14:17:00 -72.186473 0.006659
FIRELBFGS: 148 14:17:00 -72.186475 0.004181
FIRELBFGS: 149 14:17:00 -72.186491 0.007610
FIRELBFGS: 150 14:17:00 -72.186485 0.009571
FIRELBFGS: 151 14:17:00 -72.186495 0.006566
FIRELBFGS: 152 14:17:00 -72.186491 0.003034
FIRELBFGS: 153 14:17:01 -72.186489 0.002764
FIRELBFGS: 154 14:17:01 -72.186499 0.005128
FIRELBFGS: 155 14:17:01 -72.186495 0.007303
FIRELBFGS: 156 14:17:01 -72.186496 0.007636
FIRELBFGS: 157 14:17:01 -72.186509 0.005736
FIRELBFGS: 158 14:17:01 -72.186513 0.005726
FIRELBFGS: 159 14:17:01 -72.186504 0.004978
FIRELBFGS: 160 14:17:01 -72.186515 0.005859
FIRELBFGS: 161 14:17:01 -72.186520 0.006614
FIRELBFGS: 162 14:17:01 -72.186517 0.004194
FIRELBFGS: 163 14:17:01 -72.186520 0.002570
FIRELBFGS: 164 14:17:01 -72.186519 0.002108
FIRELBFGS: 165 14:17:01 -72.186520 0.001654
FIRELBFGS: 166 14:17:01 -72.186518 0.001206
FIRELBFGS: 167 14:17:01 -72.186520 0.000938
次に、LBFGSで大きく振動したり、ポテンシャルを登ってしまったAr二原子の例を見てみましょう。
このような例でも破綻することなく、頑健に収束しています。
[32]:
lennard_jones_trajectory(FIRELBFGS, 1.5)
Step Time Energy fmax
FIRELBFGS: 0 14:17:07 -0.314857 1.158029
FIRELBFGS: 1 14:17:07 -0.342894 1.264380
FIRELBFGS: 2 14:17:07 -0.410024 1.512419
FIRELBFGS: 3 14:17:07 -0.546744 1.968308
FIRELBFGS: 4 14:17:07 -0.812180 2.383891
FIRELBFGS: 5 14:17:07 -0.862167 5.585844
FIRELBFGS: 6 14:17:07 -0.966377 2.149190
FIRELBFGS: 7 14:17:07 -0.991817 0.522302
FIRELBFGS: 8 14:17:07 -0.992148 0.491275
FIRELBFGS: 9 14:17:07 -0.992732 0.429940
FIRELBFGS: 10 14:17:07 -0.993427 0.339973
FIRELBFGS: 11 14:17:07 -0.994057 0.224455
FIRELBFGS: 12 14:17:07 -0.994451 0.088426
FIRELBFGS: 13 14:17:07 -0.994489 0.060610
FIRELBFGS: 14 14:17:07 -0.994490 0.059506
FIRELBFGS: 15 14:17:07 -0.994492 0.057320
FIRELBFGS: 16 14:17:07 -0.994495 0.054094
FIRELBFGS: 17 14:17:07 -0.994499 0.049889
Distance in opt trajectory: [1.5, 1.4768394233790767, 1.42839123774426, 1.349694678243147, 1.231631957294667, 1.0658914242047832, 1.0938206432760949, 1.132495812583008, 1.1318429355798734, 1.1305759647542164, 1.1287715691339621, 1.1265422074054834, 1.1240322773166505, 1.1214118150599879, 1.1214307556785812, 1.1214682920288386, 1.1215237409768002, 1.1215960942672192]
[33]:
lennard_jones_trajectory(FIRELBFGS, 1.28)
Step Time Energy fmax
FIRELBFGS: 0 14:17:10 -0.697220 2.324552
FIRELBFGS: 1 14:17:10 -0.807702 2.387400
FIRELBFGS: 2 14:17:10 -0.987241 0.821987
FIRELBFGS: 3 14:17:10 -0.520102 13.570402
FIRELBFGS: 4 14:17:10 -0.971689 1.903593
FIRELBFGS: 5 14:17:10 -0.939109 1.840004
FIRELBFGS: 6 14:17:10 -0.943288 1.793714
FIRELBFGS: 7 14:17:10 -0.951214 1.694594
FIRELBFGS: 8 14:17:10 -0.961964 1.529352
FIRELBFGS: 9 14:17:10 -0.974030 1.277952
FIRELBFGS: 10 14:17:10 -0.985244 0.914936
FIRELBFGS: 11 14:17:10 -0.992867 0.414203
FIRELBFGS: 12 14:17:10 -0.994043 0.239611
FIRELBFGS: 13 14:17:10 -0.994061 0.235001
FIRELBFGS: 14 14:17:10 -0.994095 0.225890
FIRELBFGS: 15 14:17:10 -0.994143 0.212487
FIRELBFGS: 16 14:17:10 -0.994201 0.195100
FIRELBFGS: 17 14:17:10 -0.994265 0.174121
FIRELBFGS: 18 14:17:10 -0.994330 0.150013
FIRELBFGS: 19 14:17:10 -0.994391 0.123296
FIRELBFGS: 20 14:17:10 -0.994449 0.091418
FIRELBFGS: 21 14:17:10 -0.994495 0.054286
FIRELBFGS: 22 14:17:10 -0.994519 0.012221
Distance in opt trajectory: [1.28, 1.2335089629222111, 1.1392699328939198, 1.0285911631178963, 1.0964431714345662, 1.1738131464936081, 1.1715131413850235, 1.1669709941864936, 1.160310604717144, 1.1517385248429906, 1.1415690055780352, 1.130255815977825, 1.118424872251103, 1.1184997505580012, 1.1186480666965142, 1.1188669733107373, 1.1191522821112356, 1.1194985597723826, 1.1198992503739216, 1.1203468201911215, 1.1208857684008715, 1.121520438052656, 1.1222486280893154]
更にきちんと学びたい方のために¶
以下書籍にニュートン法ファミリーや準ニュートン法、LineSearchなどの手法に関する詳細な解説があります。
“Numerical Optimization” Jorge Nocedal, Stephan Wright https://doi.org/10.1007/978-0-387-40065-5
FIREは元論文を読むと理解しやすいでしょう。
また、どの手法もASEの実装を読むことで勉強することもできます。
[34]:
from ase.optimize import FIRE
??FIRE
Init signature:
FIRE(
atoms: ase.atoms.Atoms,
restart: Optional[str] = None,
logfile: Union[IO, str] = '-',
trajectory: Optional[str] = None,
dt: float = 0.1,
maxstep: Optional[float] = None,
maxmove: Optional[float] = None,
dtmax: float = 1.0,
Nmin: int = 5,
finc: float = 1.1,
fdec: float = 0.5,
astart: float = 0.1,
fa: float = 0.99,
a: float = 0.1,
master: Optional[bool] = None,
downhill_check: bool = False,
position_reset_callback: Optional[Callable] = None,
force_consistent=<object object at 0x7ff6a0fa36e0>,
)
Docstring: Base-class for all structure optimization classes.
Source:
class FIRE(Optimizer):
@deprecated(
"Use of `maxmove` is deprecated. Use `maxstep` instead.",
category=FutureWarning,
callback=_forbid_maxmove,
)
def __init__(
self,
atoms: Atoms,
restart: Optional[str] = None,
logfile: Union[IO, str] = '-',
trajectory: Optional[str] = None,
dt: float = 0.1,
maxstep: Optional[float] = None,
maxmove: Optional[float] = None,
dtmax: float = 1.0,
Nmin: int = 5,
finc: float = 1.1,
fdec: float = 0.5,
astart: float = 0.1,
fa: float = 0.99,
a: float = 0.1,
master: Optional[bool] = None,
downhill_check: bool = False,
position_reset_callback: Optional[Callable] = None,
force_consistent=Optimizer._deprecated,
):
"""Parameters:
atoms: Atoms object
The Atoms object to relax.
restart: string
Pickle file used to store hessian matrix. If set, file with
such a name will be searched and hessian matrix stored will
be used, if the file exists.
logfile: file object or str
If *logfile* is a string, a file with that name will be opened.
Use '-' for stdout.
trajectory: string
Pickle file used to store trajectory of atomic movement.
dt: float
Initial time step. Defualt value is 0.1
maxstep: float
Used to set the maximum distance an atom can move per
iteration (default value is 0.2).
dtmax: float
Maximum time step. Default value is 1.0
Nmin: int
Number of steps to wait after the last time the dot product of
the velocity and force is negative (P in The FIRE article) before
increasing the time step. Default value is 5.
finc: float
Factor to increase the time step. Default value is 1.1
fdec: float
Factor to decrease the time step. Default value is 0.5
astart: float
Initial value of the parameter a. a is the Coefficient for
mixing the velocity and the force. Called alpha in the FIRE article.
Default value 0.1.
fa: float
Factor to decrease the parameter alpha. Default value is 0.99
a: float
Coefficient for mixing the velocity and the force. Called
alpha in the FIRE article. Default value 0.1.
master: boolean
Defaults to None, which causes only rank 0 to save files. If
set to true, this rank will save files.
downhill_check: boolean
Downhill check directly compares potential energies of subsequent
steps of the FIRE algorithm rather than relying on the current
product v*f that is positive if the FIRE dynamics moves downhill.
This can detect numerical issues where at large time steps the step
is uphill in energy even though locally v*f is positive, i.e. the
algorithm jumps over a valley because of a too large time step.
position_reset_callback: function(atoms, r, e, e_last)
Function that takes current *atoms* object, an array of position
*r* that the optimizer will revert to, current energy *e* and
energy of last step *e_last*. This is only called if e > e_last.
force_consistent: boolean or None
Use force-consistent energy calls (as opposed to the energy
extrapolated to 0 K). If force_consistent=None, uses
force-consistent energies if available in the calculator, but
falls back to force_consistent=False if not.
.. deprecated:: 3.19.3
Use of ``maxmove`` is deprecated; please use ``maxstep``.
"""
Optimizer.__init__(self, atoms, restart, logfile, trajectory,
master, force_consistent=force_consistent)
self.dt = dt
self.Nsteps = 0
if maxstep is not None:
self.maxstep = maxstep
else:
self.maxstep = self.defaults["maxstep"]
self.dtmax = dtmax
self.Nmin = Nmin
self.finc = finc
self.fdec = fdec
self.astart = astart
self.fa = fa
self.a = a
self.downhill_check = downhill_check
self.position_reset_callback = position_reset_callback
def initialize(self):
self.v = None
def read(self):
self.v, self.dt = self.load()
def step(self, f=None):
optimizable = self.optimizable
if f is None:
f = optimizable.get_forces()
if self.v is None:
self.v = np.zeros((len(optimizable), 3))
if self.downhill_check:
self.e_last = optimizable.get_potential_energy()
self.r_last = optimizable.get_positions().copy()
self.v_last = self.v.copy()
else:
is_uphill = False
if self.downhill_check:
e = optimizable.get_potential_energy()
# Check if the energy actually decreased
if e > self.e_last:
# If not, reset to old positions...
if self.position_reset_callback is not None:
self.position_reset_callback(
optimizable, self.r_last, e,
self.e_last)
optimizable.set_positions(self.r_last)
is_uphill = True
self.e_last = optimizable.get_potential_energy()
self.r_last = optimizable.get_positions().copy()
self.v_last = self.v.copy()
vf = np.vdot(f, self.v)
if vf > 0.0 and not is_uphill:
self.v = (1.0 - self.a) * self.v + self.a * f / np.sqrt(
np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v))
if self.Nsteps > self.Nmin:
self.dt = min(self.dt * self.finc, self.dtmax)
self.a *= self.fa
self.Nsteps += 1
else:
self.v[:] *= 0.0
self.a = self.astart
self.dt *= self.fdec
self.Nsteps = 0
self.v += self.dt * f
dr = self.dt * self.v
normdr = np.sqrt(np.vdot(dr, dr))
if normdr > self.maxstep:
dr = self.maxstep * dr / normdr
r = optimizable.get_positions()
optimizable.set_positions(r + dr)
self.dump((self.v, self.dt))
File: ~/.py39/lib/python3.9/site-packages/ase/optimize/fire.py
Type: type
Subclasses:
[ ]: