We will use the quasi-harmonic approximation (via PyIron's implementation of the popular phonopy package) to evaluate look at thermal expansion and self-diffusion in Aluminium
# Generic imports
from pyiron import Project
import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
import seaborn as sns
pr = Project("PhonopyExample")
pot = '2009--Mendelev-M-I--Al-Mg--LAMMPS--ipr1'
Because repeating code is evil.
def make_phonopy_job(template_job, name):
"""
Create a phonopy job from a reference job.
Args:
template_job (pyiron job): The job to copy.
name (str): What to call this new job.
Returns:
A new phonopy job.
"""
project = template_job.project
# What I want:
# job_type = template_job.job_type
# What I have to do instead:
job_type = pr.job_type.Lammps
ref = project.create_job(job_type, name + "_ref")
ref.structure = template_job.get_final_structure().copy()
ref.potential = template_job.potential
phono = project.create_job(pr.job_type.PhonopyJob, name)
phono.ref_job = ref
return phono
def scale_array(arr, scaler=None, new_range=1.):
"""
Linearly transforms an array so that values equal to the minimum and maximum of the
`scaler` array are mapped to the range (0, `new_range`). Note that rescaled values can
still lie outside this range if the orignal values of `arr` are outside the bounds of
`scaler`.
Args:
arr (np.array): Array to rescale.
scaler (np.array): Array by which to rescale. Default is `arr`.
new_range (float): New value for data which was the size `np.amax(scaler)`.
Default is 1.
"""
if scaler is None:
scaler = arr
return new_range * (arr - np.amin(scaler)) / np.ptp(scaler)
What does the QHA say the lattice constant is as a function of temperature?
pr_te = pr.create_group("ThermalExpansion")
If we were doing VASP instead it would be important to do the least computation as possible, so here we'll start by relaxing a simple unit cell to turn into a supercell later.
job_unit = pr_te.create_job(pr.job_type.Lammps, "UnitCell")
basis = pr_te.create_structure("Al", "fcc", 4.04)
job_unit.structure = basis
job_unit.potential = pot
job_unit.calc_minimize(pressure=0.0)
job_unit.run()
The job UnitCell was saved and received the ID: 8886147
basis_rel = job_unit.get_final_structure()
A relaxation which should take zero steps given our starting position!
job_bulk_1 = pr_te.create_job(pr.job_type.Lammps, "Bulk_1")
# The _1 here refers to the fact that the volume has been rescaled by a factor of "1.0"
# (i.e. it hasn't been rescaled)
n_reps = 3
job_bulk_1.structure = basis_rel.repeat(rep=n_reps)
job_bulk_1.potential = pot
job_bulk_1.structure.plot3d();
job_bulk_1.calc_minimize(pressure=0.0)
job_bulk_1.run()
The job Bulk_1 was saved and received the ID: 8886148
Run phonopy on the bulk supercell
phono_bulk_1 = make_phonopy_job(job_bulk_1, "PhonoBulk_1")
phono_bulk_1.run()
# Run performs a whole bunch of child calculations
# Each one has the positions slightly deformed in the symmetry-appropriate ways needed
# to get the phonon properties
The job PhonoBulk_1 was saved and received the ID: 8886149 The job supercell_phonon_0 was saved and received the ID: 8886150
# Let's see what we got...
T_min = 0
T_max = 800 # a bit below melting
T_step = 25
temperatures = np.linspace(T_min, T_max, int((T_max - T_min) / T_step))
tp_bulk_1 = phono_bulk_1.get_thermal_properties(temperatures=temperatures)
# `get_thermal_properties` uses the displacements and forces to generate phonon information
U_bulk_1 = job_bulk_1.output.energy_pot[-1]
Fvib_bulk_1 = tp_bulk_1.free_energies
plt.plot(temperatures, U_bulk_1 + Fvib_bulk_1)
plt.xlabel("Temperature [K]")
plt.ylabel("Free energy ($U+F_{vib}$) [eV]")
Text(0, 0.5, 'Free energy ($U+F_{vib}$) [eV]')
Above we have the (QHA approximation to the) free energy as a function of temperature at a fixed volume. To evaluate the thermal expansivity, we need to create the entire F(V,T) surface. To get this, we just loop over jobs like the above, but scaled to have different lattice constants.
# According to Wikipedia, the thermal expansivity is about 0.0023% / Kelvin
# So at our maximum temperature, we expect around 1.8% expansion
scale_min = -0.005
scale_max = 0.02
scale_step = 0.002
scales = np.linspace(scale_min, scale_max, int((scale_max - scale_min) / scale_step))
# Let's keep things clean by making another sub-directory
pr_scales = pr_te.create_group("ScanScales")
# Loop the phonon calculation over all the volumes
sc_bulk_rel = job_bulk_1.get_final_structure()
bulk_free_energies = np.zeros((len(scales), len(temperatures)))
for i, scale in enumerate(scales):
name_tail = "_{}".format(str(scale).replace(".", "c").replace('-', 'm'))
# Make a bulk job with the rescaled structure
# (already relaxed, by symmetry won't change, calc static will be enough)
job_bulk = pr_scales.create_job(pr.job_type.Lammps, "Bulk" + name_tail)
job_bulk.potential = pot
job_bulk.structure = sc_bulk_rel.apply_strain(epsilon=scale, return_box=True)
job_bulk.calc_static()
job_bulk.run()
U = job_bulk.output.energy_tot[-1]
# Use that job as a reference for a phonopy job
phono_bulk = make_phonopy_job(job_bulk, "PhonoBulk" + name_tail)
phono_bulk.run()
tp_bulk = phono_bulk.get_thermal_properties(temperatures=temperatures)
Fvib = tp_bulk.free_energies
# Fill in the row of free energies for this volume
bulk_free_energies[i] = Fvib + U
The job Bulk_m0c005 was saved and received the ID: 8886151 The job PhonoBulk_m0c005 was saved and received the ID: 8886152 The job supercell_phonon_0 was saved and received the ID: 8886153 The job Bulk_m0c002727272727272727 was saved and received the ID: 8886154 The job PhonoBulk_m0c002727272727272727 was saved and received the ID: 8886155 The job supercell_phonon_0 was saved and received the ID: 8886156 The job Bulk_m0c000454545454545454 was saved and received the ID: 8886157 The job PhonoBulk_m0c000454545454545454 was saved and received the ID: 8886158 The job supercell_phonon_0 was saved and received the ID: 8886159 The job Bulk_0c0018181818181818186 was saved and received the ID: 8886160 The job PhonoBulk_0c0018181818181818186 was saved and received the ID: 8886161 The job supercell_phonon_0 was saved and received the ID: 8886162 The job Bulk_0c004090909090909092 was saved and received the ID: 8886163 The job PhonoBulk_0c004090909090909092 was saved and received the ID: 8886164 The job supercell_phonon_0 was saved and received the ID: 8886165 The job Bulk_0c006363636363636366 was saved and received the ID: 8886166 The job PhonoBulk_0c006363636363636366 was saved and received the ID: 8886167 The job supercell_phonon_0 was saved and received the ID: 8886168 The job Bulk_0c008636363636363636 was saved and received the ID: 8886169 The job PhonoBulk_0c008636363636363636 was saved and received the ID: 8886170 The job supercell_phonon_0 was saved and received the ID: 8886171 The job Bulk_0c01090909090909091 was saved and received the ID: 8886172 The job PhonoBulk_0c01090909090909091 was saved and received the ID: 8886173 The job supercell_phonon_0 was saved and received the ID: 8886174 The job Bulk_0c013181818181818183 was saved and received the ID: 8886175 The job PhonoBulk_0c013181818181818183 was saved and received the ID: 8886176 The job supercell_phonon_0 was saved and received the ID: 8886177 The job Bulk_0c015454545454545457 was saved and received the ID: 8886178 The job PhonoBulk_0c015454545454545457 was saved and received the ID: 8886179 The job supercell_phonon_0 was saved and received the ID: 8886180 The job Bulk_0c01772727272727273 was saved and received the ID: 8886182 The job PhonoBulk_0c01772727272727273 was saved and received the ID: 8886183 The job supercell_phonon_0 was saved and received the ID: 8886184 The job Bulk_0c02 was saved and received the ID: 8886186 The job PhonoBulk_0c02 was saved and received the ID: 8886187 The job supercell_phonon_0 was saved and received the ID: 8886188
# The lattice constant is probably a more informative value than the 0K-relative strain
latts = basis_rel.cell[0][0] * scales
# At each temperature, find the optimal volume by a simple quadratic fit
# ...Wait, which order fit will be good enough? Let's just spot-check
free_en = bulk_free_energies[:, -1]
plt.plot(latts, free_en, color='b', label='data')
# We'll plot the fit on a much denser mesh
fit_deg = 4
p = np.polyfit(latts, free_en, deg=fit_deg)
dense_latts = np.linspace(np.amin(latts), np.amax(latts), 1000)
#dense_latts = np.linspace(0, 10, 1000)
plt.plot(dense_latts, np.polyval(p=p, x=dense_latts), color='r', label='fit')
plt.xlabel('Lattice constant [$\mathrm{\AA}$]')
plt.ylabel('Bulk free energy [eV/supercell]')
plt.legend()
# Ok, a fourth-order fit seems perfectly reasonable
<matplotlib.legend.Legend at 0x2b151bd2ac10>
# Now find optimal temperatures
best_latts = np.zeros(len(temperatures))
best_latt_guess = basis_rel.cell[0][0]
for i, T in enumerate(temperatures):
free_en = bulk_free_energies[:, i]
p = np.polyfit(latts, free_en, deg=fit_deg)
extrema = np.roots(np.polyder(p, m=1)).real # Find where first-derivative is zero
best_latts[i] = extrema[np.argmin(np.abs(extrema - best_latt_guess))]
# Check that they're resonable
print(best_latt_guess, '\n', best_latts)
4.045270475668763 [0.11555233 0.11352406 0.10694882 0.10163624 0.09885196 1.43573459 0.77014253 0.60322527 0.51918649 0.46683335 0.43047109 0.40346556 0.38247104 0.36559688 0.35168556 0.33998477 0.32998232 0.32131628 0.31372298 0.30700538 0.30101301 0.29562879 0.29076015 0.28633285 0.28228657 0.2785718 0.27514744 0.27197909 0.2690377 0.26629857 0.26374056 0.26134543]
# Let's look at the landscape
fig, ax = plt.subplots()
sns.heatmap(bulk_free_energies, ax=ax, cmap="coolwarm",
xticklabels=['{:,.0f}'.format(T) for T in temperatures],
yticklabels=['{:,.2f}'.format(a) for a in latts])
ax.set_xlabel("Temperature [K]")
ax.set_ylabel("Lattice constant [$\mathrm{\AA}$]")
# Overlaying the optimal path takes a couple changes of variables
# since the heatmap is plotting integer cells
ax.plot(scale_array(temperatures, new_range=len(temperatures)),
scale_array(best_latts, scaler=latts, new_range=len(latts)),
color='k')
[<matplotlib.lines.Line2D at 0x2b1519223a90>]
Another common use of QHA is to calculate the pre-factor for migration in a diffusion event.
In particular, the diffusion jump barrier looks like $\omega_0 = \nu_0^\star \exp(-H_\mathrm{m} / k_\mathrm{B} T)$, where $\nu_0^\star = \prod_{i=1}^{3N-3} \nu_i^\mathrm{IS} / \prod_{i=1}^{3N-4} \nu_i^\mathrm{TS}$, with IS and TS indicating the initial and transition states, respectively. Note that the transition state is missing a single frequency, which is from the instability of the transition state. It's either an imaginary mode, which I think means a negative frequency. Meanwhile, $H_\mathrm{m}$ is the enthalpic barrier (difference between the initial and transition states) and $k_\mathrm{B} T$ is the usual thermal energy term.
Typically, these sorts of investigations use the nudged elastic band (NEB) to find the 0K transition state. You can do that with our new flexible jobs, but we'll save that for later. For now we'll just "approximate" the transition state with a simple linear interpolation.
Let's start by generating and relaxing the initial and final states
pr_vac = pr.create_group("Vacancies")
# Find two adjacent sites
print(job_bulk_1.structure.positions[0])
print(job_bulk_1.structure.positions[1])
# Yep, 1 and 2 will do
[0. 0. 0.] [ 2.02263524e+00 2.02263524e+00 -7.63052415e-33]
job_vac_i = pr_vac.create_job(pr.job_type.Lammps, "VacancyInitial")
job_vac_f = pr_vac.create_job(pr.job_type.Lammps, "VacancyFinal")
job_vac_i.potential = pot
job_vac_f.potential = pot
sc_vac_i = sc_bulk_rel.copy()
sc_vac_i.pop(0)
job_vac_i.structure = sc_vac_i
sc_vac_f = sc_bulk_rel.copy()
sc_vac_f.pop(1)
job_vac_f.structure = sc_vac_f
# Relax the new vacancy structures
job_vac_i.calc_minimize(pressure=0.0)
job_vac_i.run()
job_vac_f.calc_minimize(pressure=0.0)
job_vac_f.run()
The job VacancyInitial was saved and received the ID: 8886189 The job VacancyFinal was saved and received the ID: 8886190
The pyiron implementation of phonopy makes it very easy to look at the DOS. Let's see what the effect is of introducing a vacancy, and confirm that our two vacancies are equivalent.
phon_vac_i = make_phonopy_job(job_vac_i, "PhonoVacInitial")
phon_vac_f = make_phonopy_job(job_vac_f, "PhonoVacFinal")
phon_vac_i.run()
tp_vac_i = phon_vac_i.get_thermal_properties(temperatures=temperatures)
phon_vac_f.run()
tp_vac_f = phon_vac_i.get_thermal_properties(temperatures=temperatures)
# Note that the vacancy structures spawn many more child processes
# This is because the vacancy structure has lower symmetry
The job PhonoVacInitial was saved and received the ID: 8886191 The job supercell_phonon_0 was saved and received the ID: 8886192 The job supercell_phonon_1 was saved and received the ID: 8886193 The job supercell_phonon_2 was saved and received the ID: 8886194 The job supercell_phonon_3 was saved and received the ID: 8886195 The job supercell_phonon_4 was saved and received the ID: 8886196 The job supercell_phonon_5 was saved and received the ID: 8886197 The job supercell_phonon_6 was saved and received the ID: 8886198 The job supercell_phonon_7 was saved and received the ID: 8886199 The job supercell_phonon_8 was saved and received the ID: 8886200 The job supercell_phonon_9 was saved and received the ID: 8886201 The job supercell_phonon_10 was saved and received the ID: 8886202 The job supercell_phonon_11 was saved and received the ID: 8886203 The job supercell_phonon_12 was saved and received the ID: 8886204 The job supercell_phonon_13 was saved and received the ID: 8886205 The job supercell_phonon_14 was saved and received the ID: 8886207 The job supercell_phonon_15 was saved and received the ID: 8886208 The job supercell_phonon_16 was saved and received the ID: 8886209 The job supercell_phonon_17 was saved and received the ID: 8886210 The job supercell_phonon_18 was saved and received the ID: 8886211 The job supercell_phonon_19 was saved and received the ID: 8886212 The job supercell_phonon_20 was saved and received the ID: 8886213 The job PhonoVacFinal was saved and received the ID: 8886214 The job supercell_phonon_0 was saved and received the ID: 8886215 The job supercell_phonon_1 was saved and received the ID: 8886216 The job supercell_phonon_2 was saved and received the ID: 8886217 The job supercell_phonon_3 was saved and received the ID: 8886218 The job supercell_phonon_4 was saved and received the ID: 8886219 The job supercell_phonon_5 was saved and received the ID: 8886220 The job supercell_phonon_6 was saved and received the ID: 8886221 The job supercell_phonon_7 was saved and received the ID: 8886222 The job supercell_phonon_8 was saved and received the ID: 8886223 The job supercell_phonon_9 was saved and received the ID: 8886225 The job supercell_phonon_10 was saved and received the ID: 8886226 The job supercell_phonon_11 was saved and received the ID: 8886227 The job supercell_phonon_12 was saved and received the ID: 8886228 The job supercell_phonon_13 was saved and received the ID: 8886229 The job supercell_phonon_14 was saved and received the ID: 8886230 The job supercell_phonon_15 was saved and received the ID: 8886231 The job supercell_phonon_16 was saved and received the ID: 8886232 The job supercell_phonon_17 was saved and received the ID: 8886233 The job supercell_phonon_18 was saved and received the ID: 8886234 The job supercell_phonon_19 was saved and received the ID: 8886235 The job supercell_phonon_20 was saved and received the ID: 8886236
fig, ax = plt.subplots()
phono_bulk_1.plot_dos(ax=ax, color='b', label='bulk')
phon_vac_i.plot_dos(ax=ax, color='r', label='vac_i')
phon_vac_f.plot_dos(ax=ax, color='orange', label='vac_f')
plt.legend()
<matplotlib.legend.Legend at 0x2b157cd72890>
Now we get the attempt frequency by comparing the individual phonon spectra of initial and transition states
# Interpolate initial and final positions to guesstimate the transition state
sc_vac_ts = sc_vac_i.copy()
sc_vac_ts.positions = 0.5 * (sc_vac_i.positions + sc_vac_f.positions)
job_vac_ts = pr_vac.create_job(pr.job_type.Lammps, "VacancyTransition")
job_vac_ts.potential = pot
job_vac_ts.structure = sc_vac_ts
# We _don't_ relax this job, or it would fall into the initial or final state!
job_vac_ts.calc_static()
job_vac_ts.run()
The job VacancyTransition was saved and received the ID: 8886237
phon_vac_ts = make_phonopy_job(job_vac_ts, "PhonoVacTransition")
phon_vac_ts.run()
tp_vac_ts = phon_vac_ts.get_thermal_properties(temperatures=temperatures)
The job PhonoVacTransition was saved and received the ID: 8886238 The job supercell_phonon_0 was saved and received the ID: 8886239 The job supercell_phonon_1 was saved and received the ID: 8886240 The job supercell_phonon_2 was saved and received the ID: 8886241 The job supercell_phonon_3 was saved and received the ID: 8886242 The job supercell_phonon_4 was saved and received the ID: 8886243 The job supercell_phonon_5 was saved and received the ID: 8886244 The job supercell_phonon_6 was saved and received the ID: 8886245 The job supercell_phonon_7 was saved and received the ID: 8886246 The job supercell_phonon_8 was saved and received the ID: 8886247 The job supercell_phonon_9 was saved and received the ID: 8886248 The job supercell_phonon_10 was saved and received the ID: 8886249 The job supercell_phonon_11 was saved and received the ID: 8886250 The job supercell_phonon_12 was saved and received the ID: 8886251 The job supercell_phonon_13 was saved and received the ID: 8886252 The job supercell_phonon_14 was saved and received the ID: 8886253 The job supercell_phonon_15 was saved and received the ID: 8886254 The job supercell_phonon_16 was saved and received the ID: 8886255 The job supercell_phonon_17 was saved and received the ID: 8886256 The job supercell_phonon_18 was saved and received the ID: 8886257 The job supercell_phonon_19 was saved and received the ID: 8886258 The job supercell_phonon_20 was saved and received the ID: 8886259 The job supercell_phonon_21 was saved and received the ID: 8886260 The job supercell_phonon_22 was saved and received the ID: 8886261 The job supercell_phonon_23 was saved and received the ID: 8886262 The job supercell_phonon_24 was saved and received the ID: 8886263 The job supercell_phonon_25 was saved and received the ID: 8886264 The job supercell_phonon_26 was saved and received the ID: 8886265 The job supercell_phonon_27 was saved and received the ID: 8886267 The job supercell_phonon_28 was saved and received the ID: 8886269 The job supercell_phonon_29 was saved and received the ID: 8886270 The job supercell_phonon_30 was saved and received the ID: 8886271 The job supercell_phonon_31 was saved and received the ID: 8886272 The job supercell_phonon_32 was saved and received the ID: 8886273 The job supercell_phonon_33 was saved and received the ID: 8886274 The job supercell_phonon_34 was saved and received the ID: 8886275 The job supercell_phonon_35 was saved and received the ID: 8886276 The job supercell_phonon_36 was saved and received the ID: 8886277 The job supercell_phonon_37 was saved and received the ID: 8886278 The job supercell_phonon_38 was saved and received the ID: 8886279 The job supercell_phonon_39 was saved and received the ID: 8886280 The job supercell_phonon_40 was saved and received the ID: 8886281 The job supercell_phonon_41 was saved and received the ID: 8886282 The job supercell_phonon_42 was saved and received the ID: 8886283 The job supercell_phonon_43 was saved and received the ID: 8886284 The job supercell_phonon_44 was saved and received the ID: 8886285 The job supercell_phonon_45 was saved and received the ID: 8886286 The job supercell_phonon_46 was saved and received the ID: 8886287 The job supercell_phonon_47 was saved and received the ID: 8886288 The job supercell_phonon_48 was saved and received the ID: 8886289 The job supercell_phonon_49 was saved and received the ID: 8886290 The job supercell_phonon_50 was saved and received the ID: 8886291 The job supercell_phonon_51 was saved and received the ID: 8886292 The job supercell_phonon_52 was saved and received the ID: 8886293 The job supercell_phonon_53 was saved and received the ID: 8886294 The job supercell_phonon_54 was saved and received the ID: 8886295 The job supercell_phonon_55 was saved and received the ID: 8886296 The job supercell_phonon_56 was saved and received the ID: 8886297 The job supercell_phonon_57 was saved and received the ID: 8886298 The job supercell_phonon_58 was saved and received the ID: 8886299 The job supercell_phonon_59 was saved and received the ID: 8886300 The job supercell_phonon_60 was saved and received the ID: 8886301 The job supercell_phonon_61 was saved and received the ID: 8886302 The job supercell_phonon_62 was saved and received the ID: 8886303 The job supercell_phonon_63 was saved and received the ID: 8886304 The job supercell_phonon_64 was saved and received the ID: 8886305 The job supercell_phonon_65 was saved and received the ID: 8886306 The job supercell_phonon_66 was saved and received the ID: 8886307 The job supercell_phonon_67 was saved and received the ID: 8886308 The job supercell_phonon_68 was saved and received the ID: 8886309 The job supercell_phonon_69 was saved and received the ID: 8886310 The job supercell_phonon_70 was saved and received the ID: 8886311 The job supercell_phonon_71 was saved and received the ID: 8886312 The job supercell_phonon_72 was saved and received the ID: 8886313 The job supercell_phonon_73 was saved and received the ID: 8886314 The job supercell_phonon_74 was saved and received the ID: 8886315 The job supercell_phonon_75 was saved and received the ID: 8886316 The job supercell_phonon_76 was saved and received the ID: 8886317 The job supercell_phonon_77 was saved and received the ID: 8886318 The job supercell_phonon_78 was saved and received the ID: 8886319 The job supercell_phonon_79 was saved and received the ID: 8886320 The job supercell_phonon_80 was saved and received the ID: 8886321 The job supercell_phonon_81 was saved and received the ID: 8886322 The job supercell_phonon_82 was saved and received the ID: 8886323 The job supercell_phonon_83 was saved and received the ID: 8886324 The job supercell_phonon_84 was saved and received the ID: 8886325 The job supercell_phonon_85 was saved and received the ID: 8886326 The job supercell_phonon_86 was saved and received the ID: 8886327 The job supercell_phonon_87 was saved and received the ID: 8886328 The job supercell_phonon_88 was saved and received the ID: 8886329 The job supercell_phonon_89 was saved and received the ID: 8886330 The job supercell_phonon_90 was saved and received the ID: 8886331 The job supercell_phonon_91 was saved and received the ID: 8886332 The job supercell_phonon_92 was saved and received the ID: 8886333 The job supercell_phonon_93 was saved and received the ID: 8886334 The job supercell_phonon_94 was saved and received the ID: 8886335
# The transition state has an imaginary mode (frequency < 0), let's see it
fig, ax = plt.subplots()
phon_vac_i.plot_dos(ax=ax, color='r', label='initial')
phon_vac_ts.plot_dos(ax=ax, color='b', label='transition')
plt.legend()
<matplotlib.legend.Legend at 0x2b157c09cc90>
To calculate the attempt frequency, we'll ignore both the negative mode of the transition state (which we were warned about in the equation), as well as the three frequencies which correspond to rigid translation and are very near zero, and sometimes dip to be negative. Phonopy sorts the frequencies by magnitude, so we can just skip the first three and four for the initial and transition states, respectively. We take them at q=0.
freq_i = phon_vac_i.phonopy.get_frequencies(0)[3:]
freq_ts = phon_vac_i.phonopy.get_frequencies(0)[4:]
print(np.prod(freq_i))
6.870293244293476e+236
Recall: $\nu_0^\star = \prod_{i=1}^{3N-3} \nu_i^\mathrm{IS} / \prod_{i=1}^{3N-4} \nu_i^\mathrm{TS}$
# Products are dangerous beasts, so we'll do a little numeric magic
nu = np.prod(freq_i[:-1] / freq_ts) * freq_i[-1]
print("Attempt frequency is ", nu, "THz (10^-12 s)")
Attempt frequency is 2.6826762430167848 THz (10^-12 s)
Mantina et al. (PRL 2008) report $\nu = 19.3$ THz using DFT and NEB, so our linearly-interpolated "transition state" with EAM is actually not doing so poorly.
There are many more things you can do with phonopy, including looking directly at the force constants, the Hessian matrix, etc. But hopefully this is a useful starting point.