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.18.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.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)
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
[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()
Let’s run a structural optimization of this simple potential form by starting the distance away from the stability point, 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()
The first case is when FIRE is used. The optimum value of 1.12A is reached in a stable manner.
[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]
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.
[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]
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.
[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]
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.
[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]
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.
[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 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.
[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]
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!
[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 |
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.)
[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"])
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.
[28]:
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.
[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
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.
[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]
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.
[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:
[ ]: