Phonopy in pyiron

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_atomistics 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"

Helper functions

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_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.0):
    """
    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)

Thermal Expansion

What does the QHA say the lattice constant is as a function of temperature?

pr_te = pr.create_group("ThermalExpansion")

Relax the unit cell

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)

/var/folders/9p/rztyv06d0xv4h26cyv8nrw3m0000gq/T/ipykernel_60014/2777024177.py:1: DeprecationWarning: pyiron_atomistics.project.create_structure is deprecated: Use create.structure.* methods instead.
  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: 2

basis_rel = job_unit.get_structure()

Relax the bulk supercell

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: 3

Calculate phonons

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: 4
The job PhonoBulk_1_ref_0 was saved and received the ID: 5

# 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]')

Calculate thermal expansivity

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.002
scale_max = 0.025
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_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_m0c002 was saved and received the ID: 6
The job PhonoBulk_m0c002 was saved and received the ID: 7
The job PhonoBulk_m0c002_ref_0 was saved and received the ID: 8
The job Bulk_0c0002500000000000002 was saved and received the ID: 9
The job PhonoBulk_0c0002500000000000002 was saved and received the ID: 10
The job PhonoBulk_0c0002500000000000002_ref_0 was saved and received the ID: 11
The job Bulk_0c0025000000000000005 was saved and received the ID: 12
The job PhonoBulk_0c0025000000000000005 was saved and received the ID: 13
The job PhonoBulk_0c0025000000000000005_ref_0 was saved and received the ID: 14
The job Bulk_0c004750000000000001 was saved and received the ID: 15
The job PhonoBulk_0c004750000000000001 was saved and received the ID: 16
The job PhonoBulk_0c004750000000000001_ref_0 was saved and received the ID: 17
The job Bulk_0c007000000000000001 was saved and received the ID: 18
The job PhonoBulk_0c007000000000000001 was saved and received the ID: 19
The job PhonoBulk_0c007000000000000001_ref_0 was saved and received the ID: 20
The job Bulk_0c009250000000000001 was saved and received the ID: 21
The job PhonoBulk_0c009250000000000001 was saved and received the ID: 22
The job PhonoBulk_0c009250000000000001_ref_0 was saved and received the ID: 23
The job Bulk_0c011500000000000002 was saved and received the ID: 24
The job PhonoBulk_0c011500000000000002 was saved and received the ID: 25
The job PhonoBulk_0c011500000000000002_ref_0 was saved and received the ID: 26
The job Bulk_0c01375 was saved and received the ID: 27
The job PhonoBulk_0c01375 was saved and received the ID: 28
The job PhonoBulk_0c01375_ref_0 was saved and received the ID: 29
The job Bulk_0c016 was saved and received the ID: 30
The job PhonoBulk_0c016 was saved and received the ID: 31
The job PhonoBulk_0c016_ref_0 was saved and received the ID: 32
The job Bulk_0c018250000000000002 was saved and received the ID: 33
The job PhonoBulk_0c018250000000000002 was saved and received the ID: 34
The job PhonoBulk_0c018250000000000002_ref_0 was saved and received the ID: 35
The job Bulk_0c020500000000000004 was saved and received the ID: 36
The job PhonoBulk_0c020500000000000004 was saved and received the ID: 37
The job PhonoBulk_0c020500000000000004_ref_0 was saved and received the ID: 38
The job Bulk_0c02275 was saved and received the ID: 39
The job PhonoBulk_0c02275 was saved and received the ID: 40
The job PhonoBulk_0c02275_ref_0 was saved and received the ID: 41
The job Bulk_0c025 was saved and received the ID: 42
The job PhonoBulk_0c025 was saved and received the ID: 43
The job PhonoBulk_0c025_ref_0 was saved and received the ID: 44

# The lattice constant is probably a more informative value than the 0K-relative strain
latts = basis_rel.cell[0][0] * (1 + 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 0x16da333d0>

# 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 
 [4.05959821 4.05963134 4.0600255  4.06101957 4.06247856 4.0642089
 4.06606872 4.06796789 4.0698518  4.07168845 4.07345997 4.07515736
 4.07677712 4.07831912 4.0797853  4.08117875 4.0825032  4.08376262
 4.08496107 4.08610249 4.08719068 4.08822921 4.08922145 4.09017051
 4.09107929 4.09195045 4.09278647 4.09358962 4.09436199 4.09510551
 4.09582196 4.09651298]

# 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 0x155720090>]

Vacancies and diffusion

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.

Stable vacancy structures

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: 45
The job VacancyFinal was saved and received the ID: 46

DOS

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: 47
The job PhonoVacInitial_ref_0 was saved and received the ID: 48
The job PhonoVacInitial_ref_1 was saved and received the ID: 49
The job PhonoVacInitial_ref_2 was saved and received the ID: 50
The job PhonoVacInitial_ref_3 was saved and received the ID: 51
The job PhonoVacInitial_ref_4 was saved and received the ID: 52
The job PhonoVacInitial_ref_5 was saved and received the ID: 53
The job PhonoVacInitial_ref_6 was saved and received the ID: 54
The job PhonoVacInitial_ref_7 was saved and received the ID: 55
The job PhonoVacInitial_ref_8 was saved and received the ID: 56
The job PhonoVacInitial_ref_9 was saved and received the ID: 57
The job PhonoVacInitial_ref_10 was saved and received the ID: 58
The job PhonoVacInitial_ref_11 was saved and received the ID: 59
The job PhonoVacInitial_ref_12 was saved and received the ID: 60
The job PhonoVacInitial_ref_13 was saved and received the ID: 61
The job PhonoVacInitial_ref_14 was saved and received the ID: 62
The job PhonoVacInitial_ref_15 was saved and received the ID: 63
The job PhonoVacInitial_ref_16 was saved and received the ID: 64
The job PhonoVacInitial_ref_17 was saved and received the ID: 65
The job PhonoVacInitial_ref_18 was saved and received the ID: 66
The job PhonoVacInitial_ref_19 was saved and received the ID: 67
The job PhonoVacInitial_ref_20 was saved and received the ID: 68
The job PhonoVacFinal was saved and received the ID: 69
The job PhonoVacFinal_ref_0 was saved and received the ID: 70
The job PhonoVacFinal_ref_1 was saved and received the ID: 71
The job PhonoVacFinal_ref_2 was saved and received the ID: 72
The job PhonoVacFinal_ref_3 was saved and received the ID: 73
The job PhonoVacFinal_ref_4 was saved and received the ID: 74
The job PhonoVacFinal_ref_5 was saved and received the ID: 75
The job PhonoVacFinal_ref_6 was saved and received the ID: 76
The job PhonoVacFinal_ref_7 was saved and received the ID: 77
The job PhonoVacFinal_ref_8 was saved and received the ID: 78
The job PhonoVacFinal_ref_9 was saved and received the ID: 79
The job PhonoVacFinal_ref_10 was saved and received the ID: 80
The job PhonoVacFinal_ref_11 was saved and received the ID: 81
The job PhonoVacFinal_ref_12 was saved and received the ID: 82
The job PhonoVacFinal_ref_13 was saved and received the ID: 83
The job PhonoVacFinal_ref_14 was saved and received the ID: 84
The job PhonoVacFinal_ref_15 was saved and received the ID: 85
The job PhonoVacFinal_ref_16 was saved and received the ID: 86
The job PhonoVacFinal_ref_17 was saved and received the ID: 87
The job PhonoVacFinal_ref_18 was saved and received the ID: 88
The job PhonoVacFinal_ref_19 was saved and received the ID: 89
The job PhonoVacFinal_ref_20 was saved and received the ID: 90

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()

/var/folders/9p/rztyv06d0xv4h26cyv8nrw3m0000gq/T/ipykernel_60014/3735091431.py:2: DeprecationWarning: pyiron_atomistics.atomistics.master.phonopy.plot_dos(ax=Axes(0.125,0.11;0.775x0.77)) is deprecated.
  phono_bulk_1.plot_dos(ax=ax, color='b', label='bulk')
/var/folders/9p/rztyv06d0xv4h26cyv8nrw3m0000gq/T/ipykernel_60014/3735091431.py:3: DeprecationWarning: pyiron_atomistics.atomistics.master.phonopy.plot_dos(ax=Axes(0.125,0.11;0.775x0.77)) is deprecated.
  phon_vac_i.plot_dos(ax=ax, color='r', label='vac_i')
/var/folders/9p/rztyv06d0xv4h26cyv8nrw3m0000gq/T/ipykernel_60014/3735091431.py:4: DeprecationWarning: pyiron_atomistics.atomistics.master.phonopy.plot_dos(ax=Axes(0.125,0.11;0.775x0.77)) is deprecated.
  phon_vac_f.plot_dos(ax=ax, color='orange', label='vac_f')

<matplotlib.legend.Legend at 0x16d71dc50>

Attempt frequency

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: 91

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: 92
The job PhonoVacTransition_ref_0 was saved and received the ID: 93
The job PhonoVacTransition_ref_1 was saved and received the ID: 94
The job PhonoVacTransition_ref_2 was saved and received the ID: 95
The job PhonoVacTransition_ref_3 was saved and received the ID: 96
The job PhonoVacTransition_ref_4 was saved and received the ID: 97
The job PhonoVacTransition_ref_5 was saved and received the ID: 98
The job PhonoVacTransition_ref_6 was saved and received the ID: 99
The job PhonoVacTransition_ref_7 was saved and received the ID: 100
The job PhonoVacTransition_ref_8 was saved and received the ID: 101
The job PhonoVacTransition_ref_9 was saved and received the ID: 102
The job PhonoVacTransition_ref_10 was saved and received the ID: 103
The job PhonoVacTransition_ref_11 was saved and received the ID: 104
The job PhonoVacTransition_ref_12 was saved and received the ID: 105
The job PhonoVacTransition_ref_13 was saved and received the ID: 106
The job PhonoVacTransition_ref_14 was saved and received the ID: 107
The job PhonoVacTransition_ref_15 was saved and received the ID: 108
The job PhonoVacTransition_ref_16 was saved and received the ID: 109
The job PhonoVacTransition_ref_17 was saved and received the ID: 110
The job PhonoVacTransition_ref_18 was saved and received the ID: 111
The job PhonoVacTransition_ref_19 was saved and received the ID: 112
The job PhonoVacTransition_ref_20 was saved and received the ID: 113
The job PhonoVacTransition_ref_21 was saved and received the ID: 114
The job PhonoVacTransition_ref_22 was saved and received the ID: 115
The job PhonoVacTransition_ref_23 was saved and received the ID: 116
The job PhonoVacTransition_ref_24 was saved and received the ID: 117
The job PhonoVacTransition_ref_25 was saved and received the ID: 118
The job PhonoVacTransition_ref_26 was saved and received the ID: 119
The job PhonoVacTransition_ref_27 was saved and received the ID: 120
The job PhonoVacTransition_ref_28 was saved and received the ID: 121
The job PhonoVacTransition_ref_29 was saved and received the ID: 122
The job PhonoVacTransition_ref_30 was saved and received the ID: 123
The job PhonoVacTransition_ref_31 was saved and received the ID: 124
The job PhonoVacTransition_ref_32 was saved and received the ID: 125
The job PhonoVacTransition_ref_33 was saved and received the ID: 126
The job PhonoVacTransition_ref_34 was saved and received the ID: 127
The job PhonoVacTransition_ref_35 was saved and received the ID: 128
The job PhonoVacTransition_ref_36 was saved and received the ID: 129
The job PhonoVacTransition_ref_37 was saved and received the ID: 130
The job PhonoVacTransition_ref_38 was saved and received the ID: 131
The job PhonoVacTransition_ref_39 was saved and received the ID: 132
The job PhonoVacTransition_ref_40 was saved and received the ID: 133
The job PhonoVacTransition_ref_41 was saved and received the ID: 134
The job PhonoVacTransition_ref_42 was saved and received the ID: 135
The job PhonoVacTransition_ref_43 was saved and received the ID: 136
The job PhonoVacTransition_ref_44 was saved and received the ID: 137
The job PhonoVacTransition_ref_45 was saved and received the ID: 138
The job PhonoVacTransition_ref_46 was saved and received the ID: 139
The job PhonoVacTransition_ref_47 was saved and received the ID: 140
The job PhonoVacTransition_ref_48 was saved and received the ID: 141
The job PhonoVacTransition_ref_49 was saved and received the ID: 142
The job PhonoVacTransition_ref_50 was saved and received the ID: 143
The job PhonoVacTransition_ref_51 was saved and received the ID: 144
The job PhonoVacTransition_ref_52 was saved and received the ID: 145
The job PhonoVacTransition_ref_53 was saved and received the ID: 146
The job PhonoVacTransition_ref_54 was saved and received the ID: 147
The job PhonoVacTransition_ref_55 was saved and received the ID: 148
The job PhonoVacTransition_ref_56 was saved and received the ID: 149
The job PhonoVacTransition_ref_57 was saved and received the ID: 150
The job PhonoVacTransition_ref_58 was saved and received the ID: 151
The job PhonoVacTransition_ref_59 was saved and received the ID: 152
The job PhonoVacTransition_ref_60 was saved and received the ID: 153
The job PhonoVacTransition_ref_61 was saved and received the ID: 154
The job PhonoVacTransition_ref_62 was saved and received the ID: 155
The job PhonoVacTransition_ref_63 was saved and received the ID: 156
The job PhonoVacTransition_ref_64 was saved and received the ID: 157
The job PhonoVacTransition_ref_65 was saved and received the ID: 158
The job PhonoVacTransition_ref_66 was saved and received the ID: 159
The job PhonoVacTransition_ref_67 was saved and received the ID: 160
The job PhonoVacTransition_ref_68 was saved and received the ID: 161
The job PhonoVacTransition_ref_69 was saved and received the ID: 162
The job PhonoVacTransition_ref_70 was saved and received the ID: 163
The job PhonoVacTransition_ref_71 was saved and received the ID: 164
The job PhonoVacTransition_ref_72 was saved and received the ID: 165
The job PhonoVacTransition_ref_73 was saved and received the ID: 166
The job PhonoVacTransition_ref_74 was saved and received the ID: 167
The job PhonoVacTransition_ref_75 was saved and received the ID: 168
The job PhonoVacTransition_ref_76 was saved and received the ID: 169
The job PhonoVacTransition_ref_77 was saved and received the ID: 170
The job PhonoVacTransition_ref_78 was saved and received the ID: 171
The job PhonoVacTransition_ref_79 was saved and received the ID: 172
The job PhonoVacTransition_ref_80 was saved and received the ID: 173
The job PhonoVacTransition_ref_81 was saved and received the ID: 174
The job PhonoVacTransition_ref_82 was saved and received the ID: 175
The job PhonoVacTransition_ref_83 was saved and received the ID: 176
The job PhonoVacTransition_ref_84 was saved and received the ID: 177
The job PhonoVacTransition_ref_85 was saved and received the ID: 178
The job PhonoVacTransition_ref_86 was saved and received the ID: 179
The job PhonoVacTransition_ref_87 was saved and received the ID: 180
The job PhonoVacTransition_ref_88 was saved and received the ID: 181
The job PhonoVacTransition_ref_89 was saved and received the ID: 182
The job PhonoVacTransition_ref_90 was saved and received the ID: 183
The job PhonoVacTransition_ref_91 was saved and received the ID: 184
The job PhonoVacTransition_ref_92 was saved and received the ID: 185
The job PhonoVacTransition_ref_93 was saved and received the ID: 186
The job PhonoVacTransition_ref_94 was saved and received the ID: 187

# 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()

/var/folders/9p/rztyv06d0xv4h26cyv8nrw3m0000gq/T/ipykernel_60014/2028323302.py:3: DeprecationWarning: pyiron_atomistics.atomistics.master.phonopy.plot_dos(ax=Axes(0.125,0.11;0.775x0.77)) is deprecated.
  phon_vac_i.plot_dos(ax=ax, color='r', label='initial')
/var/folders/9p/rztyv06d0xv4h26cyv8nrw3m0000gq/T/ipykernel_60014/2028323302.py:4: DeprecationWarning: pyiron_atomistics.atomistics.master.phonopy.plot_dos(ax=Axes(0.125,0.11;0.775x0.77)) is deprecated.
  phon_vac_ts.plot_dos(ax=ax, color='b', label='transition')

<matplotlib.legend.Legend at 0x155961a10>

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, 0, 0))[3:]
freq_ts = phon_vac_i.phonopy.get_frequencies((0, 0, 0))[4:]

print(np.prod(freq_i))

6.870476287488258e+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.682686767439414 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.