Optimization algorithm¶
Until the previous section, we applied structural optimization to various examples. In this section, we will learn about the local optimization algorithms used in structural optimization.
Algorithm type¶
As local optimization algorithms, ASE offers the FIRE method, which performs optimization with Molecular Dynamics (MD)-like behavior, and the BFGS, LBFGS, BFGSLineSearch=QuasiNewton, and LBFGSLineSearch methods are provided as quasi-newton methods.
Algorithm |
Group |
Description |
|---|---|---|
MDMin |
MD-like |
Performs MD, but resets momentum to 0 if the inner product of momentum and force is negative |
FIRE |
MD-like |
MD, with various additional innovations to make it fast and robust. |
BFGS |
Quasi-Newton |
Approximate the Hessian from the optimization trajectory and run the Newton method using the approximated Hessian |
LBFGS |
Quasi-Newton |
BFGS method to run at high speed and low memory |
BFGSLineSearch |
Quasi-Newton |
BFGS determines the direction of optimization direction, and LineSearch determines the step size. |
LBFGSLineSearch |
Quasi-Newton |
LBFGS determines the direction of optimization direction, and LineSearch determines the step size. |
We exclude MDMin from the benchmark because of its large hyperparameter dependence and frequent structural optimization failures. Let’s look at other optimizers.
FIRE¶
FIRE basically performs MD, but some ingenious techniques are included for fast convergence.
Newton method¶
The Newton method uses the Hessian and gradient to determine the next step. The Newton method is a quadratic convergence method, which means that when the potential is close to a quadratic function (close to the minima), it approaches the exact solution at the speed of the square and is very fast. With the Hessian as \(H\) and the gradient as \(\vec{g}\), the next optimization step \(\vec{p}\) is expressed as follows
\(\vec{p} = - H^{-1} \vec{g}\)
If the potential were strictly quadratic, it would converge to a minimum in one step.
Quasi-Newton Method¶
Since the evaluation of the Hessian is onerous (for example, it can be obtained numerically by computing force 6N times), it is often more meaningful to use the time spent computing the Hessian to proceed with structural optimization. Therefore, using the optimization trajectory, a quasi-Newton method approximates the Hessian in a short time. Since the Hessian is not exact, it does not converge in as few steps as the Newton method but still converges fast.
BFGS¶
The BFGS method is the most standard Hessian approximation method among the quasi-Newtonian methods. Various other methods are known, such as the SR1 and DFP methods, but the BFGS method is implemented in ASE.
LBFGS¶
One drawback of the BFGS method is that it requires storing a matrix of the size of the variable’s dimension to optimize and calculate the matrix-vector product. This means \(O(N^2)\) for dimension \(N\) in space and \(O(N^2)\) in time are required, which is not suitable for optimizing multidimensional functions. Therefore, the LBFGS method is a method that can reproduce BFGS with a computational complexity of \(O(N)\) by using only the information from the last few trajectories and transforming it so that the results are almost the same as the BFGS method. Since the computational complexity of PFP is \(O(N)\), the computational time in Optimizer becomes noticeable from about 300 atoms using the \(O(N^2)\) method like BFGS. However, in LBFGS, the computational time in optimizer is not noticeable even if the number of atoms increases.
LineSearch¶
If the potential is not quadratic, the Newton method may oscillate. Also, the quasi-Newton method may not provide an appropriate step size because the Hessian may not be accurate. In these cases, LineSearch can be used to stabilize the optimization in order to converge the optimization. After determining the direction of optimization using the quasi-Newton method, it is desirable to decide a step size such that the energy is sufficiently small in that direction. LineSearch determines such a step size by actually calculating several points. With BFGS, there may be cases where the energy increases as a result of calculation according to the Newton method, but BFGSLineSearch searches for a good step size.
Benchmark¶
[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
Comparison of number of steps (force calculation)¶
Below is a simple example of a structural optimization of a CH3CHO molecule. The structure obtained by molecule("CH3CHO") is not strictly stable, but it is a fairly good initial structure. Therefore, when starting from such a structure, the quasi-Newton method can be used to complete the structural optimization in a much smaller number of PFP calls.
First, let’s benchmark how many times force is calculated until the structural optimization is complete. Here, we take the structure of the bulk of Pt moved randomly as the initial structure and measure the number of force calculations until the optimal structure is achieved.
[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 had the highest number of force evaluation, while BFGS and LBFGS had about the same number. The method with LineSearch also had the same number of force evaluations as BFGS and LBFGS.
Comparison of speed per step¶
Next, let us compare the speed per step of each of the optimization methods.
[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)
In BFGS and BFGSLineSearch, the calculation time per step increases as the number of atoms increases.
The BFGS method implemented in ASE performs eigenvalue decomposition at each step to find the Hessian inverse matrix. Since the Hessian has a size of 3N, where N is the number of atoms, and the eigenvalue decomposition takes \(O(N^3)\) computational time, the BFGS method in ASE slows down rapidly on the order of the cube of the number of atoms.
BFGSLineSearch directly approximates the Hessian inverse matrix without performing Hessian eigenvalue decomposition. This makes it faster than BFGS, \(O(N^2)\) computational time because the matrix-vector product becomes the rate-limiting factor.
However, it is much slower than LBFGS and FIRE, which operate on \(O(N)\) when the number of atoms is large.
Use of different optimization algorithms¶
We have seen that FIRE requires more force calls than other methods, and that BFGS and BFGSLineSearch are more rate-limiting in the internal routines of the optimizer than in force evaluation when the number of atoms is large. So in what cases should we use LBFGS and LBFGSLineSearch? We will now look at some cases where even those algorithms do not work.
Examples of LBFGS and BFGS fail¶
Let’s look at the next example. This example uses a very simple potential, the Lennard-Jones (LJ) potential, and performs structural optimization with initial values away from its stability point.
LJ potential is expressed as follows, and can be visualized as below
[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()
Let’s run a structural optimization of this simple potential form by starting the distance away from the stability point, 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()
The first case is when FIRE is used. The optimum value of 1.12A is reached in a stable manner.
[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]
If we start the BFGS method at 1.5A, we will find that it jumps too far at certain steps, resulting in unstable behavior that jumps far past the stability point. (where it seems to jump from 1.18A to 0.78A.) Thus, the distance will oscillate significantly before reaching the optimum value of 1.12A, even in this simple example. When the function surface is very different from the second-order function and the second-order approximation method is used, its estimation of minima is far from the actual minima.
[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]
The same thing happens with the LBFGS method. If such oscillations occur in a polyatomic system, the optimization may not be completed for a long time.
[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]
There are also other cases where the initial value does not reach the stability point correctly if it starts far from the stability point. This may be caused by a problem with the initialization of Hessian or other reasons.
[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]
Thus, even in cases where BFGS and LBFGS do not work, LineSearch can be used to converge correctly. LBFGSLineSearch is a better optimizer in most cases.
[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 appears to finish in a few steps, but the actual number of calculations is a bit higher because the part of LineSearch being done internally is not counted as a step. BFGSLineSearch shows Step[FC], where FC is the actual number of force evaluations.
[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]
We have seen here that there are examples where BFGS and LBFGS do not work. But even in these instances, LBFGSLineSearch has worked well.
Our conclusion so far is that LBFGSLineSearch is the best optimizer in many cases. Below we discuss the shortcomings of the other optimizers.
FIRE: Long number of steps to convergence
BFGS: Takes a long time when the number of atoms is large.
BFGSLineSearch: Takes a long time when the number of atoms is large.
LBFGS: May oscillate or increase energy sometimes.
However, it is not always the case that using LBFGSLineSearch always works, unfortunately.
Example of LBFGSLineSearch fails¶
LBFGSLineSearch is a complete method, and it will work well in many cases. However, there are rare cases where LBFGSLineSearch does not work well. For example, this is often the case when one wants to perform precise structural optimization with a small fmax. If you want to optimize the rotation of the methyl group of toluene, you need to set fmax=0.001 eV/A to reach the minimum value for the rotation. In this example, LineSearch will not converge easily, and one step will take a
long time, or the energy will not drop, or in some cases, LineSearch will fail, and an error will occur with 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 |
In this example, LBFGS converged, but FIRE did not.
Also, LBFGSLineSearch took a very long time per step in the latter half and did not converge, in some cases giving the error RuntimeError: LineSearch failed!. In most cases, LBFGSLineSearch has difficulty converging around fmax=0.002 eV/A. In some cases, even fmax=0.05 may not be reached if the shape of the potential is worse.
If LBFGSLineSearch stops working like this, you may have to re-run LBFGSLineSearch to get it to work. (Sometimes it doesn’t work.)
[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"])
Implementation of optimization algorithms based on the above findings¶
LBFGSLineSearch is a good optimizer, but it can be very slow with small fmax, and in some cases, it can also produce errors. FIRE is often slower than LBFGSLineSearch, but it is robust in many cases. However, it may not converge easily when fmax is small. LBFGS converges even with small fmax, but in some cases, it may oscillate or increase its energy, as seen in the Lennard-Jones example.
Heuristically, we want to use the second-order Newton method with good convergence close to the stability point where the quadratic function approximation is a good approximation. And we adopt LBFGS since LBFGSLineSearch is not numerically stable in some cases. For unstable structures where the quadratic approximation is not a good approximation, we adopt FIRE or LBFGSLineSearch to drop energy robustly.
However, it is difficult to predict at what fmax cases LBFGSLineSearch will not be numerically stable. FIRE is more stable when the method is applied uniformly to unknown materials, such as in high-throughput calculations.
To summarize our findings so far, the optimization algorithms that seem to work well heuristically are
In unstable places, such as the early stages of structural optimization, use FIRE to move down the gradient.
Once a certain point of stability is reached (small values of force, etc.), the LBFGS method is used to attempt fast convergence.
This algorithm is implemented as ``FIRELBFGS`` optimizer in matlantis-features.
[30]:
from matlantis_features.ase_ext.optimize import FIRELBFGS
First, let’s look at an example where LBFGSLineSearch does not work: the rotation of the methyl group of toluene. FIRELBFGS works in these examples and converges robustly, slower than LBFGS but much faster than 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
Next, let’s look at an example of an Ar diatomic that has oscillated significantly or climbed a potential in LBFGS. Even in this example, the system converges robustly without breaking down.
[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]
Further reading¶
The following books provide detailed explanations of the Newton method family, quasi-Newton methods, LineSearch, and other methods.
“Numerical Optimization” Jorge Nocedal, Stephan Wright https://doi.org/10.1007/978-0-387-40065-5
The original paper is the best source to understand FIRE.
You may also study algorithms by reading ASE implementations.
[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: