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.5.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.06 sec.
Running BFGS...
Done in 0.08 sec.
Running LBFGS...
Done in 0.06 sec.
Running BFGSLineSearch...
Done in 0.06 sec.
Running LBFGSLineSearch...
Done in 0.06 sec.
----- n_atoms 200 -----
Running FIRE...
Done in 0.08 sec.
Running BFGS...
Done in 0.16 sec.
Running LBFGS...
Done in 0.08 sec.
Running BFGSLineSearch...
Done in 0.11 sec.
Running LBFGSLineSearch...
Done in 0.08 sec.
----- n_atoms 300 -----
Running FIRE...
Done in 0.12 sec.
Running BFGS...
Done in 0.32 sec.
Running LBFGS...
Done in 0.11 sec.
Running BFGSLineSearch...
Done in 0.19 sec.
Running LBFGSLineSearch...
Done in 0.11 sec.
----- n_atoms 400 -----
Running FIRE...
Done in 0.15 sec.
Running BFGS...
Done in 0.59 sec.
Running LBFGS...
Done in 0.15 sec.
Running BFGSLineSearch...
Done in 0.31 sec.
Running LBFGSLineSearch...
Done in 0.15 sec.
----- n_atoms 500 -----
Running FIRE...
Done in 0.18 sec.
Running BFGS...
Done in 0.94 sec.
Running LBFGS...
Done in 0.18 sec.
Running BFGSLineSearch...
Done in 0.40 sec.
Running LBFGSLineSearch...
Done in 0.18 sec.
----- n_atoms 600 -----
Running FIRE...
Done in 0.21 sec.
Running BFGS...
Done in 1.58 sec.
Running LBFGS...
Done in 0.21 sec.
Running BFGSLineSearch...
Done in 0.56 sec.
Running LBFGSLineSearch...
Done in 0.21 sec.
----- n_atoms 700 -----
Running FIRE...
Done in 0.25 sec.
Running BFGS...
Done in 2.45 sec.
Running LBFGS...
Done in 0.24 sec.
Running BFGSLineSearch...
Done in 0.81 sec.
Running LBFGSLineSearch...
Done in 0.24 sec.
----- n_atoms 800 -----
Running FIRE...
Done in 0.28 sec.
Running BFGS...
Done in 3.41 sec.
Running LBFGS...
Done in 0.27 sec.
Running BFGSLineSearch...
Done in 1.03 sec.
Running LBFGSLineSearch...
Done in 0.27 sec.
----- n_atoms 900 -----
Running FIRE...
Done in 0.31 sec.
Running BFGS...
Done in 4.87 sec.
Running LBFGS...
Done in 0.30 sec.
Running BFGSLineSearch...
Done in 1.22 sec.
Running LBFGSLineSearch...
Done in 0.30 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\}\]
[15]:
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 から離れたところからスタートし構造最適化をしてみましょう。
[16]:
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 にたどり着くことが出来ています。
[17]:
lennard_jones_trajectory(FIRE, 1.5)
Step Time Energy fmax
*Force-consistent energies used in optimization.
FIRE: 0 06:34:10 -0.314857* 1.1580
FIRE: 1 06:34:10 -0.342894* 1.2644
FIRE: 2 06:34:10 -0.410024* 1.5124
FIRE: 3 06:34:10 -0.546744* 1.9683
FIRE: 4 06:34:10 -0.812180* 2.3839
FIRE: 5 06:34:10 -0.862167* 5.5858
FIRE: 6 06:34:10 -0.966377* 2.1492
FIRE: 7 06:34:10 -0.991817* 0.5223
FIRE: 8 06:34:10 -0.992148* 0.4913
FIRE: 9 06:34:10 -0.992732* 0.4299
FIRE: 10 06:34:10 -0.993427* 0.3400
FIRE: 11 06:34:10 -0.994057* 0.2245
FIRE: 12 06:34:10 -0.994451* 0.0884
FIRE: 13 06:34:10 -0.994489* 0.0606
FIRE: 14 06:34:10 -0.994490* 0.0595
FIRE: 15 06:34:10 -0.994492* 0.0573
FIRE: 16 06:34:10 -0.994495* 0.0541
FIRE: 17 06:34:10 -0.994499* 0.0499
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]
[18]:
lennard_jones_trajectory(FIRE, 1.28)
Step Time Energy fmax
*Force-consistent energies used in optimization.
FIRE: 0 06:34:10 -0.697220* 2.3246
FIRE: 1 06:34:10 -0.807702* 2.3874
FIRE: 2 06:34:10 -0.987241* 0.8220
FIRE: 3 06:34:10 -0.520102* 13.5704
FIRE: 4 06:34:10 -0.971689* 1.9036
FIRE: 5 06:34:10 -0.939109* 1.8400
FIRE: 6 06:34:10 -0.943288* 1.7937
FIRE: 7 06:34:10 -0.951214* 1.6946
FIRE: 8 06:34:10 -0.961964* 1.5294
FIRE: 9 06:34:10 -0.974030* 1.2780
FIRE: 10 06:34:10 -0.985244* 0.9149
FIRE: 11 06:34:10 -0.992867* 0.4142
FIRE: 12 06:34:10 -0.994043* 0.2396
FIRE: 13 06:34:10 -0.994061* 0.2350
FIRE: 14 06:34:10 -0.994095* 0.2259
FIRE: 15 06:34:10 -0.994143* 0.2125
FIRE: 16 06:34:10 -0.994201* 0.1951
FIRE: 17 06:34:10 -0.994265* 0.1741
FIRE: 18 06:34:10 -0.994330* 0.1500
FIRE: 19 06:34:10 -0.994391* 0.1233
FIRE: 20 06:34:10 -0.994449* 0.0914
FIRE: 21 06:34:10 -0.994495* 0.0543
FIRE: 22 06:34:10 -0.994519* 0.0122
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次関数と異なる場合に、その極小値の見積もりが大きくはずれてしまうために起こります。
[19]:
lennard_jones_trajectory(BFGS, 1.5)
Step Time Energy fmax
BFGS: 0 06:34:11 -0.314857 1.1580
BFGS: 1 06:34:11 -0.355681 1.3124
BFGS: 2 06:34:11 -0.916039 2.0410
BFGS: 3 06:34:11 55.304434 974.4875
BFGS: 4 06:34:11 -0.917740 2.0288
BFGS: 5 06:34:11 -0.919418 2.0163
BFGS: 6 06:34:11 -0.747080 8.5176
BFGS: 7 06:34:11 -0.965068 1.4729
BFGS: 8 06:34:11 -0.984660 0.9397
BFGS: 9 06:34:11 -0.992324 0.5286
BFGS: 10 06:34:11 -0.994444 0.0924
BFGS: 11 06:34:11 -0.994520 0.0073
Distance in opt trajectory: [1.5, 1.4669134619701096, 1.1856710148074785, 0.7856710148074794, 1.1848349876531552, 1.1840057048694352, 1.0494144415980853, 1.1582431525834571, 1.1421985974452118, 1.1139254943723882, 1.1241042745909056, 1.122589482834374]
同じことがLBFGS法でも起こります。
このような振動が多原子系で起こった場合、本当にいつまでたっても最適化が終了しないような場合もあります。
[20]:
lennard_jones_trajectory(LBFGS, 1.28)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 06:34:11 -0.697220* 2.3246
LBFGS: 1 06:34:11 -0.854683* 2.3149
LBFGS: 2 06:34:11 33.771592* 599.7519
LBFGS: 3 06:34:11 -0.858237* 2.3057
LBFGS: 4 06:34:11 -0.861748* 2.2959
LBFGS: 5 06:34:11 17.305333* 318.4918
LBFGS: 6 06:34:11 -0.867638* 2.2781
LBFGS: 7 06:34:11 -0.873395* 2.2589
LBFGS: 8 06:34:11 5.604073* 121.1548
LBFGS: 9 06:34:11 -0.885566* 2.2113
LBFGS: 10 06:34:11 -0.897001* 2.1571
LBFGS: 11 06:34:11 0.369118* 30.7872
LBFGS: 12 06:34:11 -0.925197* 1.9707
LBFGS: 13 06:34:11 -0.947015* 1.7491
LBFGS: 14 06:34:11 -0.915603* 4.0090
LBFGS: 15 06:34:11 -0.985429* 0.9069
LBFGS: 16 06:34:11 -0.993163* 0.3771
LBFGS: 17 06:34:11 -0.994465* 0.0806
LBFGS: 18 06:34:11 -0.994520* 0.0054
Distance in opt trajectory: [1.28, 1.213584232746016, 0.813584232746016, 1.2120462614888505, 1.2105202846657441, 0.8506725019776178, 1.2079447803538974, 1.2054073889508412, 0.9079986088524981, 1.19996383178627, 1.1947303311793396, 0.986665993291516, 1.181106905537103, 1.1694092735266515, 1.077087999344094, 1.141365556759249, 1.1295077111909042, 1.1210690913427293, 1.1225559889249783]
また、他にも初期値が安定点から離れたところからスタートした場合に正しく安定点にたどり着けていないケースも有ります。 これは、Hessianの初期化がうまく行かないなどの理由で起きていると考えられます。
[21]:
lennard_jones_trajectory(LBFGS, 1.50)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 06:34:12 -0.314857* 1.1580
LBFGS: 1 06:34:12 -0.355681* 1.3124
LBFGS: 2 06:34:12 -0.129756* 0.4473
LBFGS: 3 06:34:12 -0.079409* 0.2630
LBFGS: 4 06:34:12 -0.040472* 0.1297
LBFGS: 5 06:34:12 -0.021155* 0.0689
LBFGS: 6 06:34:12 -0.009646* 0.0357
Distance in opt trajectory: [1.5, 1.4669134619701096, 1.748155909132741, 1.8935676545167726, 2.101103976550205, 2.302941039822042, 2.531922335240822]
このようにBFGSやLBFGSがうまくいかない例であっても、LineSearchを行うことで正しく収束させることが出来ます。
LBFGSLineSearchはほとんどの場合において最良のOptimizerです。
[22]:
lennard_jones_trajectory(LBFGSLineSearch, 1.5)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGSLineSearch: 0 06:34:12 -0.314857* 1.1580
LBFGSLineSearch: 1 06:34:13 -0.994517* 0.0214
Distance in opt trajectory: [1.5, 1.1220880716401869]
LBFGSLineSearchでは1ステップで終了しているように見えますが、内部でLineSearchを行っている部分がステップにカウントされていないため、実際の計算数はもう少し多いです。 BFGSLineSearchでは Step[FC]という表示がありますが、このFCが実際のフォース評価回数です。
[23]:
lennard_jones_trajectory(BFGSLineSearch, 1.5)
Step[ FC] Time Energy fmax
*Force-consistent energies used in optimization.
BFGSLineSearch: 0[ 0] 06:34:13 -0.314857* 1.1580
BFGSLineSearch: 1[ 5] 06:34:13 -0.994057* 0.2244
BFGSLineSearch: 2[ 7] 06:34:13 -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! となってエラーになってしまう場合もあります。[24]:
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
[25]:
view_ngl(atoms_0, representations=["ball+stick"], h=600, w=600)
[26]:
# 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_22968/983813568.py in <module>
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:
~/.local/lib/python3.7/site-packages/ase/optimize/optimize.py in run(self, fmax, steps)
267 if steps:
268 self.max_steps = steps
--> 269 return Dynamics.run(self)
270
271 def converged(self, forces=None):
~/.local/lib/python3.7/site-packages/ase/optimize/optimize.py in run(self)
154 *steps*."""
155
--> 156 for converged in Dynamics.irun(self):
157 pass
158 return converged
~/.local/lib/python3.7/site-packages/ase/optimize/optimize.py in irun(self)
133
134 # compute the next step
--> 135 self.step()
136 self.nsteps += 1
137
~/.local/lib/python3.7/site-packages/ase/optimize/lbfgs.py in step(self, f)
145 if self.use_line_search is True:
146 e = self.func(r)
--> 147 self.line_search(r, g, e)
148 dr = (self.alpha_k * self.p).reshape(len(self.atoms), -1)
149 else:
~/.local/lib/python3.7/site-packages/ase/optimize/lbfgs.py in line_search(self, r, g, e)
239 c2=.46, stpmax=50.)
240 if self.alpha_k is None:
--> 241 raise RuntimeError('LineSearch failed!')
242
243
RuntimeError: LineSearch failed!
[27]:
steps["index"] = "nsteps"
df = pd.DataFrame(steps)
df = df.set_index("index")
df
[27]:
| FIRE | LBFGS | LBFGSLineSearch | |
|---|---|---|---|
| index | |||
| nsteps | 200 | 80 | 627 |
この例ではLBFGSは収束しますが、FIREは収束しませんでした。
また、LBFGSLineSearchは後半に1ステップあたりにすごく長い時間がかかった上に収束せず、場合によっては
RuntimeError: LineSearch failed! というエラーが出ます。経験的に多くの場合はfmax=0.002eV/A 付近からLBFGSLineSearchが収束しにくくなりますが、
もっとポテンシャルの形状が悪い場合はfmax=0.05に到達しないケースもあります。
ちなみに、このようにLBFGSLineSearchが進まなくなってしまった場合はLBFGSLineSearchを実行しなおすことで、うまくいくようになることがあります。（うまくいかない時もあります）
[28]:
with LBFGSLineSearch(atoms) as opt:
opt.run(0.001)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGSLineSearch: 0 06:37:34 -72.184941* 0.0481
LBFGSLineSearch: 1 06:37:34 -72.184962* 0.0229
LBFGSLineSearch: 2 06:37:34 -72.184966* 0.0086
LBFGSLineSearch: 3 06:37:34 -72.184971* 0.0084
LBFGSLineSearch: 4 06:37:34 -72.184984* 0.0089
LBFGSLineSearch: 5 06:37:34 -72.184987* 0.0082
LBFGSLineSearch: 6 06:37:34 -72.184993* 0.0051
LBFGSLineSearch: 7 06:37:35 -72.184991* 0.0059
LBFGSLineSearch: 8 06:37:35 -72.185005* 0.0060
LBFGSLineSearch: 9 06:37:35 -72.185003* 0.0060
LBFGSLineSearch: 10 06:37:35 -72.185006* 0.0059
LBFGSLineSearch: 11 06:37:36 -72.185005* 0.0059
LBFGSLineSearch: 12 06:37:36 -72.185004* 0.0059
LBFGSLineSearch: 13 06:37:37 -72.185005* 0.0059
LBFGSLineSearch: 14 06:37:37 -72.185005* 0.0059
LBFGSLineSearch: 15 06:37:37 -72.185005* 0.0059
LBFGSLineSearch: 16 06:37:38 -72.185004* 0.0059
LBFGSLineSearch: 17 06:37:38 -72.185005* 0.0097
LBFGSLineSearch: 18 06:37:38 -72.185004* 0.0097
LBFGSLineSearch: 19 06:37:39 -72.185004* 0.0097
LBFGSLineSearch: 20 06:37:39 -72.185008* 0.0092
LBFGSLineSearch: 21 06:37:39 -72.185006* 0.0092
LBFGSLineSearch: 22 06:37:40 -72.185007* 0.0092
LBFGSLineSearch: 23 06:37:40 -72.185006* 0.0092
LBFGSLineSearch: 24 06:37:40 -72.185007* 0.0092
LBFGSLineSearch: 25 06:37:41 -72.185008* 0.0092
LBFGSLineSearch: 26 06:37:41 -72.185007* 0.0092
LBFGSLineSearch: 27 06:37:41 -72.185008* 0.0074
LBFGSLineSearch: 28 06:37:42 -72.185011* 0.0075
LBFGSLineSearch: 29 06:37:42 -72.185014* 0.0079
LBFGSLineSearch: 30 06:37:42 -72.185013* 0.0079
LBFGSLineSearch: 31 06:37:42 -72.185012* 0.0079
LBFGSLineSearch: 32 06:37:43 -72.185010* 0.0122
LBFGSLineSearch: 33 06:37:43 -72.185012* 0.0122
LBFGSLineSearch: 34 06:37:43 -72.185012* 0.0122
LBFGSLineSearch: 35 06:37:43 -72.185011* 0.0122
LBFGSLineSearch: 36 06:37:44 -72.185011* 0.0122
LBFGSLineSearch: 37 06:37:44 -72.185009* 0.0122
LBFGSLineSearch: 38 06:37:45 -72.185011* 0.0122
LBFGSLineSearch: 39 06:37:45 -72.185011* 0.0122
LBFGSLineSearch: 40 06:37:45 -72.185010* 0.0122
LBFGSLineSearch: 41 06:37:46 -72.185010* 0.0122
LBFGSLineSearch: 42 06:37:46 -72.185009* 0.0122
LBFGSLineSearch: 43 06:37:47 -72.185008* 0.0121
LBFGSLineSearch: 44 06:37:47 -72.185015* 0.0136
LBFGSLineSearch: 45 06:37:47 -72.185015* 0.0136
LBFGSLineSearch: 46 06:37:47 -72.185014* 0.0136
LBFGSLineSearch: 47 06:37:48 -72.185014* 0.0136
LBFGSLineSearch: 48 06:37:48 -72.185015* 0.0136
LBFGSLineSearch: 49 06:37:48 -72.185014* 0.0136
LBFGSLineSearch: 50 06:37:49 -72.185015* 0.0136
LBFGSLineSearch: 51 06:37:49 -72.185016* 0.0136
LBFGSLineSearch: 52 06:37:49 -72.185015* 0.0136
LBFGSLineSearch: 53 06:37:49 -72.185015* 0.0136
LBFGSLineSearch: 54 06:37:50 -72.185015* 0.0136
LBFGSLineSearch: 55 06:37:50 -72.185015* 0.0136
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/tmp/ipykernel_22968/1423403903.py in <module>
1 with LBFGSLineSearch(atoms) as opt:
----> 2 opt.run(0.001)
~/.local/lib/python3.7/site-packages/ase/optimize/optimize.py in run(self, fmax, steps)
267 if steps:
268 self.max_steps = steps
--> 269 return Dynamics.run(self)
270
271 def converged(self, forces=None):
~/.local/lib/python3.7/site-packages/ase/optimize/optimize.py in run(self)
154 *steps*."""
155
--> 156 for converged in Dynamics.irun(self):
157 pass
158 return converged
~/.local/lib/python3.7/site-packages/ase/optimize/optimize.py in irun(self)
133
134 # compute the next step
--> 135 self.step()
136 self.nsteps += 1
137
~/.local/lib/python3.7/site-packages/ase/optimize/lbfgs.py in step(self, f)
145 if self.use_line_search is True:
146 e = self.func(r)
--> 147 self.line_search(r, g, e)
148 dr = (self.alpha_k * self.p).reshape(len(self.atoms), -1)
149 else:
~/.local/lib/python3.7/site-packages/ase/optimize/lbfgs.py in line_search(self, r, g, e)
239 c2=.46, stpmax=50.)
240 if self.alpha_k is None:
--> 241 raise RuntimeError('LineSearch failed!')
242
243
RuntimeError: LineSearch failed!
[29]:
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 として実装をしています。
[30]:
from matlantis_features.ase_ext.optimize import FIRELBFGS
まずはLBFGSLineSearchが破綻した例であるトルエンのメチル基の回転を見てみましょう。
このような例でも破綻することなく、LBFGSよりは遅いがFIREよりはずっと速く、頑健に収束しています。
[31]:
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
[32]:
view_ngl(atoms_0, representations=["ball+stick"], h=600, w=600)
[33]:
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 06:38:11 -72.042262 1.9298
FIRELBFGS: 1 06:38:11 -72.115648 0.6318
FIRELBFGS: 2 06:38:11 -72.123152 1.4667
FIRELBFGS: 3 06:38:11 -72.136884 1.1383
FIRELBFGS: 4 06:38:11 -72.153114 0.5753
FIRELBFGS: 5 06:38:11 -72.160545 0.3966
FIRELBFGS: 6 06:38:11 -72.159824 0.6265
FIRELBFGS: 7 06:38:11 -72.160903 0.5933
FIRELBFGS: 8 06:38:11 -72.162829 0.5287
FIRELBFGS: 9 06:38:11 -72.165230 0.4361
FIRELBFGS: 10 06:38:11 -72.167638 0.3211
FIRELBFGS: 11 06:38:11 -72.169681 0.2252
FIRELBFGS: 12 06:38:11 -72.171167 0.2068
FIRELBFGS: 13 06:38:11 -72.172165 0.2283
FIRELBFGS: 14 06:38:11 -72.173021 0.2531
FIRELBFGS: 15 06:38:12 -72.174056 0.3027
FIRELBFGS: 16 06:38:12 -72.175576 0.3459
FIRELBFGS: 17 06:38:12 -72.177591 0.3275
FIRELBFGS: 18 06:38:12 -72.179676 0.2362
FIRELBFGS: 19 06:38:12 -72.181151 0.1671
FIRELBFGS: 20 06:38:12 -72.181706 0.1748
FIRELBFGS: 21 06:38:12 -72.181787 0.1623
FIRELBFGS: 22 06:38:12 -72.181925 0.1381
FIRELBFGS: 23 06:38:12 -72.182081 0.1042
FIRELBFGS: 24 06:38:12 -72.182244 0.0692
FIRELBFGS: 25 06:38:12 -72.182364 0.0661
FIRELBFGS: 26 06:38:12 -72.182445 0.0625
FIRELBFGS: 27 06:38:12 -72.182485 0.0689
FIRELBFGS: 28 06:38:12 -72.182522 0.0930
FIRELBFGS: 29 06:38:12 -72.182577 0.1027
FIRELBFGS: 30 06:38:12 -72.182690 0.0954
FIRELBFGS: 31 06:38:12 -72.182858 0.0733
FIRELBFGS: 32 06:38:12 -72.183047 0.0498
FIRELBFGS: 33 06:38:12 -72.183208 0.0403
FIRELBFGS: 34 06:38:12 -72.183328 0.0712
FIRELBFGS: 35 06:38:13 -72.183456 0.0721
FIRELBFGS: 36 06:38:13 -72.183613 0.0480
FIRELBFGS: 37 06:38:13 -72.183761 0.0392
FIRELBFGS: 38 06:38:13 -72.183848 0.0549
FIRELBFGS: 39 06:38:13 -72.183912 0.0463
FIRELBFGS: 40 06:38:13 -72.183962 0.0496
FIRELBFGS: 41 06:38:13 -72.183982 0.0402
FIRELBFGS: 42 06:38:13 -72.184009 0.0269
FIRELBFGS: 43 06:38:13 -72.184026 0.0257
FIRELBFGS: 44 06:38:13 -72.184041 0.0248
FIRELBFGS: 45 06:38:13 -72.184050 0.0312
FIRELBFGS: 46 06:38:13 -72.184071 0.0308
FIRELBFGS: 47 06:38:13 -72.184102 0.0226
FIRELBFGS: 48 06:38:13 -72.184136 0.0205
FIRELBFGS: 49 06:38:13 -72.184159 0.0176
FIRELBFGS: 50 06:38:13 -72.184184 0.0151
FIRELBFGS: 51 06:38:13 -72.184221 0.0193
FIRELBFGS: 52 06:38:13 -72.184407 0.0294
FIRELBFGS: 53 06:38:13 -72.184473 0.0249
FIRELBFGS: 54 06:38:13 -72.184519 0.0271
FIRELBFGS: 55 06:38:14 -72.184575 0.0225
FIRELBFGS: 56 06:38:14 -72.184658 0.0358
FIRELBFGS: 57 06:38:14 -72.184724 0.0253
FIRELBFGS: 58 06:38:14 -72.184755 0.0163
FIRELBFGS: 59 06:38:14 -72.184776 0.0111
FIRELBFGS: 60 06:38:14 -72.184779 0.0073
FIRELBFGS: 61 06:38:14 -72.184780 0.0070
FIRELBFGS: 62 06:38:14 -72.184790 0.0074
FIRELBFGS: 63 06:38:14 -72.184807 0.0108
FIRELBFGS: 64 06:38:14 -72.184813 0.0098
FIRELBFGS: 65 06:38:14 -72.184829 0.0080
FIRELBFGS: 66 06:38:14 -72.184830 0.0085
FIRELBFGS: 67 06:38:14 -72.184842 0.0101
FIRELBFGS: 68 06:38:14 -72.184852 0.0098
FIRELBFGS: 69 06:38:14 -72.184858 0.0057
FIRELBFGS: 70 06:38:14 -72.184859 0.0037
FIRELBFGS: 71 06:38:14 -72.184866 0.0037
FIRELBFGS: 72 06:38:14 -72.184866 0.0044
FIRELBFGS: 73 06:38:14 -72.184865 0.0041
FIRELBFGS: 74 06:38:15 -72.184870 0.0037
FIRELBFGS: 75 06:38:15 -72.184871 0.0033
FIRELBFGS: 76 06:38:15 -72.184870 0.0050
FIRELBFGS: 77 06:38:15 -72.184879 0.0053
FIRELBFGS: 78 06:38:15 -72.184881 0.0047
FIRELBFGS: 79 06:38:15 -72.184877 0.0041
FIRELBFGS: 80 06:38:15 -72.184878 0.0052
FIRELBFGS: 81 06:38:15 -72.184881 0.0050
FIRELBFGS: 82 06:38:15 -72.184882 0.0043
FIRELBFGS: 83 06:38:15 -72.184887 0.0043
FIRELBFGS: 84 06:38:15 -72.184886 0.0058
FIRELBFGS: 85 06:38:15 -72.184897 0.0094
FIRELBFGS: 86 06:38:15 -72.184904 0.0128
FIRELBFGS: 87 06:38:15 -72.184920 0.0128
FIRELBFGS: 88 06:38:15 -72.184937 0.0134
FIRELBFGS: 89 06:38:15 -72.184962 0.0138
FIRELBFGS: 90 06:38:15 -72.184996 0.0207
FIRELBFGS: 91 06:38:16 -72.185054 0.0323
FIRELBFGS: 92 06:38:16 -72.185137 0.0384
FIRELBFGS: 93 06:38:16 -72.185245 0.0380
FIRELBFGS: 94 06:38:16 -72.185407 0.0407
FIRELBFGS: 95 06:38:16 -72.185615 0.0433
FIRELBFGS: 96 06:38:16 -72.185314 0.1413
FIRELBFGS: 97 06:38:16 -72.185818 0.0851
FIRELBFGS: 98 06:38:16 -72.186076 0.0979
FIRELBFGS: 99 06:38:16 -72.186119 0.0839
FIRELBFGS: 100 06:38:16 -72.186170 0.0584
FIRELBFGS: 101 06:38:16 -72.186209 0.0263
FIRELBFGS: 102 06:38:16 -72.186233 0.0245
FIRELBFGS: 103 06:38:16 -72.186247 0.0331
FIRELBFGS: 104 06:38:16 -72.186263 0.0447
FIRELBFGS: 105 06:38:16 -72.186278 0.0435
FIRELBFGS: 106 06:38:16 -72.186292 0.0298
FIRELBFGS: 107 06:38:16 -72.186297 0.0245
FIRELBFGS: 108 06:38:16 -72.186296 0.0225
FIRELBFGS: 109 06:38:16 -72.186300 0.0187
FIRELBFGS: 110 06:38:17 -72.186310 0.0135
FIRELBFGS: 111 06:38:17 -72.186307 0.0106
FIRELBFGS: 112 06:38:17 -72.186311 0.0092
FIRELBFGS: 113 06:38:17 -72.186327 0.0053
FIRELBFGS: 114 06:38:17 -72.186331 0.0054
FIRELBFGS: 115 06:38:17 -72.186338 0.0074
FIRELBFGS: 116 06:38:17 -72.186339 0.0071
FIRELBFGS: 117 06:38:17 -72.186348 0.0053
FIRELBFGS: 118 06:38:17 -72.186348 0.0036
FIRELBFGS: 119 06:38:17 -72.186351 0.0038
FIRELBFGS: 120 06:38:17 -72.186353 0.0041
FIRELBFGS: 121 06:38:17 -72.186362 0.0045
FIRELBFGS: 122 06:38:17 -72.186352 0.0056
FIRELBFGS: 123 06:38:17 -72.186363 0.0038
FIRELBFGS: 124 06:38:17 -72.186363 0.0050
FIRELBFGS: 125 06:38:17 -72.186365 0.0053
FIRELBFGS: 126 06:38:17 -72.186371 0.0042
FIRELBFGS: 127 06:38:17 -72.186369 0.0047
FIRELBFGS: 128 06:38:17 -72.186375 0.0044
FIRELBFGS: 129 06:38:18 -72.186380 0.0041
FIRELBFGS: 130 06:38:18 -72.186370 0.0046
FIRELBFGS: 131 06:38:18 -72.186387 0.0057
FIRELBFGS: 132 06:38:18 -72.186383 0.0064
FIRELBFGS: 133 06:38:18 -72.186396 0.0054
FIRELBFGS: 134 06:38:18 -72.186402 0.0066
FIRELBFGS: 135 06:38:18 -72.186406 0.0087
FIRELBFGS: 136 06:38:18 -72.186414 0.0096
FIRELBFGS: 137 06:38:18 -72.186424 0.0104
FIRELBFGS: 138 06:38:18 -72.186449 0.0099
FIRELBFGS: 139 06:38:18 -72.186461 0.0083
FIRELBFGS: 140 06:38:18 -72.186464 0.0089
FIRELBFGS: 141 06:38:18 -72.186471 0.0093
FIRELBFGS: 142 06:38:18 -72.186482 0.0088
FIRELBFGS: 143 06:38:18 -72.186496 0.0101
FIRELBFGS: 144 06:38:18 -72.186506 0.0083
FIRELBFGS: 145 06:38:18 -72.186509 0.0041
FIRELBFGS: 146 06:38:18 -72.186510 0.0043
FIRELBFGS: 147 06:38:18 -72.186517 0.0042
FIRELBFGS: 148 06:38:19 -72.186515 0.0055
FIRELBFGS: 149 06:38:19 -72.186522 0.0058
FIRELBFGS: 150 06:38:19 -72.186518 0.0052
FIRELBFGS: 151 06:38:19 -72.186530 0.0057
FIRELBFGS: 152 06:38:19 -72.186532 0.0059
FIRELBFGS: 153 06:38:19 -72.186538 0.0060
FIRELBFGS: 154 06:38:19 -72.186540 0.0056
FIRELBFGS: 155 06:38:19 -72.186550 0.0069
FIRELBFGS: 156 06:38:19 -72.186546 0.0061
FIRELBFGS: 157 06:38:19 -72.186549 0.0037
FIRELBFGS: 158 06:38:19 -72.186553 0.0026
FIRELBFGS: 159 06:38:19 -72.186552 0.0022
FIRELBFGS: 160 06:38:19 -72.186551 0.0017
FIRELBFGS: 161 06:38:19 -72.186555 0.0012
FIRELBFGS: 162 06:38:19 -72.186556 0.0010
次に、LBFGSで大きく振動したり、ポテンシャルを登ってしまったAr二原子の例を見てみましょう。
このような例でも破綻することなく、頑健に収束しています。
[34]:
lennard_jones_trajectory(FIRELBFGS, 1.5)
Step Time Energy fmax
FIRELBFGS: 0 06:38:19 -0.314857 1.1580
FIRELBFGS: 1 06:38:19 -0.342894 1.2644
FIRELBFGS: 2 06:38:19 -0.410024 1.5124
FIRELBFGS: 3 06:38:19 -0.546744 1.9683
FIRELBFGS: 4 06:38:19 -0.812180 2.3839
FIRELBFGS: 5 06:38:19 -0.862167 5.5858
FIRELBFGS: 6 06:38:19 -0.966377 2.1492
FIRELBFGS: 7 06:38:19 -0.991817 0.5223
FIRELBFGS: 8 06:38:19 -0.992148 0.4913
FIRELBFGS: 9 06:38:19 -0.992732 0.4299
FIRELBFGS: 10 06:38:19 -0.993427 0.3400
FIRELBFGS: 11 06:38:19 -0.994057 0.2245
FIRELBFGS: 12 06:38:19 -0.994451 0.0884
FIRELBFGS: 13 06:38:19 -0.994489 0.0606
FIRELBFGS: 14 06:38:19 -0.994490 0.0595
FIRELBFGS: 15 06:38:19 -0.994492 0.0573
FIRELBFGS: 16 06:38:19 -0.994495 0.0541
FIRELBFGS: 17 06:38:19 -0.994499 0.0499
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]
[35]:
lennard_jones_trajectory(FIRELBFGS, 1.28)
Step Time Energy fmax
FIRELBFGS: 0 06:38:20 -0.697220 2.3246
FIRELBFGS: 1 06:38:20 -0.807702 2.3874
FIRELBFGS: 2 06:38:20 -0.987241 0.8220
FIRELBFGS: 3 06:38:20 -0.520102 13.5704
FIRELBFGS: 4 06:38:20 -0.971689 1.9036
FIRELBFGS: 5 06:38:20 -0.939109 1.8400
FIRELBFGS: 6 06:38:20 -0.943288 1.7937
FIRELBFGS: 7 06:38:20 -0.951214 1.6946
FIRELBFGS: 8 06:38:20 -0.961964 1.5294
FIRELBFGS: 9 06:38:20 -0.974030 1.2780
FIRELBFGS: 10 06:38:20 -0.985244 0.9149
FIRELBFGS: 11 06:38:20 -0.992867 0.4142
FIRELBFGS: 12 06:38:20 -0.994043 0.2396
FIRELBFGS: 13 06:38:20 -0.994061 0.2350
FIRELBFGS: 14 06:38:20 -0.994095 0.2259
FIRELBFGS: 15 06:38:20 -0.994143 0.2125
FIRELBFGS: 16 06:38:20 -0.994201 0.1951
FIRELBFGS: 17 06:38:20 -0.994265 0.1741
FIRELBFGS: 18 06:38:20 -0.994330 0.1500
FIRELBFGS: 19 06:38:20 -0.994391 0.1233
FIRELBFGS: 20 06:38:20 -0.994449 0.0914
FIRELBFGS: 21 06:38:20 -0.994495 0.0543
FIRELBFGS: 22 06:38:20 -0.994519 0.0122
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の実装を読むことで勉強することもできます。
[36]:
from ase.optimize import FIRE
??FIRE
Init signature:
FIRE(
atoms,
restart=None,
logfile='-',
trajectory=None,
dt=0.1,
maxstep=None,
maxmove=None,
dtmax=1.0,
Nmin=5,
finc=1.1,
fdec=0.5,
astart=0.1,
fa=0.99,
a=0.1,
master=None,
downhill_check=False,
position_reset_callback=None,
force_consistent=None,
)
Docstring: Base-class for all structure optimization classes.
Source:
class FIRE(Optimizer):
def __init__(self, atoms, restart=None, logfile='-', trajectory=None,
dt=0.1, maxstep=None, maxmove=None, dtmax=1.0, Nmin=5,
finc=1.1, fdec=0.5,
astart=0.1, fa=0.99, a=0.1, master=None, downhill_check=False,
position_reset_callback=None, force_consistent=None):
"""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.
trajectory: string
Pickle file used to store trajectory of atomic movement.
logfile: file object or str
If *logfile* is a string, a file with that name will be opened.
Use '-' for stdout.
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). By default (force_consistent=None) uses
force-consistent energies if available in the calculator, but
falls back to force_consistent=False if not. Only meaningful
when downhill_check is True.
"""
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
elif maxmove is not None:
self.maxstep = maxmove
warnings.warn('maxmove is deprecated; please use maxstep',
np.VisibleDeprecationWarning)
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):
atoms = self.atoms
if f is None:
f = atoms.get_forces()
if self.v is None:
self.v = np.zeros((len(atoms), 3))
if self.downhill_check:
self.e_last = atoms.get_potential_energy(
force_consistent=self.force_consistent)
self.r_last = atoms.get_positions().copy()
self.v_last = self.v.copy()
else:
is_uphill = False
if self.downhill_check:
e = atoms.get_potential_energy(
force_consistent=self.force_consistent)
# 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(atoms, self.r_last, e,
self.e_last)
atoms.set_positions(self.r_last)
is_uphill = True
self.e_last = atoms.get_potential_energy(
force_consistent=self.force_consistent)
self.r_last = atoms.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 = atoms.get_positions()
atoms.set_positions(r + dr)
self.dump((self.v, self.dt))
File: ~/.local/lib/python3.7/site-packages/ase/optimize/fire.py
Type: type
Subclasses: