AlLi phase diagram #
In this notebook, we will use the EAM potential fitted in the previous day and move towards the AlLi phase diagram. We will use many of the methods and tools we discussed in the last sessions and put them together for calculation of phase diagrams. We start with the phase diagram of AlLi from Ref. [1].
In the last session, we calculated the melting temperatures of the pure phases Al and Li, thereby arriving at two points on the phase diagram. In this notebook, we will start with the left side of the phase diagram, until \(X_{Li} < 0.5\).
As always, we start by importing the necessary modules.
import matplotlib.pyplot as plt
import numpy as np
from pyiron_atomistics import Project
from helpers import *
from calphy.integrators import kb
from tqdm.notebook import tqdm
/home/menon/miniconda3/envs/pyiron-atom-dev/lib/python3.10/site-packages/paramiko-2.10.3-py3.10.egg/paramiko/transport.py:236: CryptographyDeprecationWarning: Blowfish has been deprecated
"class": algorithms.Blowfish,
Most of the interesting features of the phase diagram at composition of \(X_{Li} < 0.5\) lies between the temperature range of 800-1000 K. Therefore, we will calculate the free energy in this range. Similar to the previous session, we will use reversible scaling to obtain the free energy in this temperature range in a single calculation. We will recalculate the free energy of pure Al in FCC lattice and pure Al liquid first.
pr = Project("phase_diagram")
Pure Al #
Solid #
We start by creating an FCC structure. The converged lattice constant has been provided so as to speed up the calculations.
structure = pr.create.structure.ase.bulk('Al', cubic=True, a=4.135).repeat(4)
We can visualise the structure.
structure.plot3d()
Our FCC lattice consists of a 5x5x5 supercell consisting of 500 atoms. Now we create a Calphy job in pyiron and assign a potential.
job_sol = pr.create.job.Calphy("xp_sol")
job_sol.structure = structure
The potential has been configured for pure Al.
job_sol.potential = "Al-atomicrex"
We have assigned the potential and structure. To speed up the calculations, we will run it on 2 cores.
job_sol.server.cores = 2
Now let’s use the calc_free_energy method.
job_sol.calc_free_energy(temperature=[800, 1000],
pressure=0,
reference_phase="solid",
n_equilibration_steps=5000,
n_switching_steps=5000)
Before we actually run the calculation, let us discuss the various parameters. temperature keyword gives the temperature range over which the free energy is to be calculated. Since we provide [700, 1000], the free energy is calculated between this range. pressure denotes the pressure of the calculation, we chose 0 in this case. Since we are using a solid FCC lattice, we set reference_phase to "solid". This means that the Einstein crystal will be used as the reference system. Finally, we have n_equilibration_steps and n_switching_steps. n_equilibration_steps denotes the number of MD steps over which the system is equilibrated to the required temperature and pressure. n_switching_steps are the number of MD steps over which the system is continuously transformed between the given interatomic potential, and the reference Einstein crystal.
Finally we run the calculation.
job_sol.run()
2022-06-07 20:54:40,399 - pyiron_log - WARNING - The job xp_sol is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
Liquid #
Before we look at the output of the previous calculation, we will also calculate the free energy of the liquid phase. For this we can use the same structure as the solid. The Calphy workflow will first superheat the structure, melt it, and then equilibrate to the required temperature and pressure. Therefore the input for the pyiron job looks fairly same.
job_lqd = pr.create.job.Calphy("xp_lqd")
job_lqd.structure = structure
job_lqd.potential = "Al-atomicrex"
job_lqd.server.cores = 2
job_lqd.calc_free_energy(temperature=[800, 1000],
pressure=0,
reference_phase="liquid",
n_equilibration_steps=5000,
n_switching_steps=5000)
The major change in the input is that the reference_phase is "liquid", instead of "solid". In this case, the Uhlenbeck-Ford model is used as the reference system instead of the Einstein crystal. Now run the job,
job_lqd.run()
2022-06-07 20:54:42,201 - pyiron_log - WARNING - The job xp_lqd is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
Now we can look at the output; and plot the free energies of the two phases and calculate the melting temperature.
pr
{'groups': ['lial_thermodynamics_composition', 'phase_diagram'], 'nodes': ['xp_alli', 'xp_sol', 'xp_lqd', 'x0_sol', 'x0_lqd', 'x0_alli', 'x1_sol', 'x1_lqd', 'x1_alli', 'x2_sol', 'x2_lqd', 'x2_alli', 'x3_sol', 'x3_lqd', 'x3_alli', 'x4_sol', 'x4_lqd', 'x4_alli']}
plt.plot(pr['xp_sol'].output.temperature, pr['xp_sol'].output.energy_free,
label="Al solid", color='#C62828')
plt.plot(pr['xp_lqd'].output.temperature, pr['xp_lqd'].output.energy_free,
label="Al liquid", color='#006899')
[<matplotlib.lines.Line2D at 0x7f622020f7f0>]
Great, finally, we will also calculate the free energy of the Al Li structure. Let’s read in the structrue from the given file and plot it.
AlLi = pr.create.structure.ase.read('LiAl_poscar', format='vasp')
AlLi.plot3d()
We are calculating the free energy at zero percent Li, there will replace all the Li atoms with Al.
AlLi[:] = 'Al'
AlLi.plot3d()
Now we run the calculation to calculate the free energy at the same temperature range.
job_AlLi = pr.create.job.Calphy("xp_alli")
job_AlLi.structure = AlLi
job_AlLi.potential = "Al-atomicrex"
job_AlLi.server.cores = 2
job_AlLi.calc_free_energy(temperature=[800, 1000],
pressure=0,
reference_phase="solid",
n_equilibration_steps=5000,
n_switching_steps=5000)
job_AlLi.run()
2022-06-07 20:54:48,944 - pyiron_log - WARNING - The job xp_alli is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
Free energy with composition #
Now we will calculate the free energy of FCC solid, liquid and B32 phases with increasing Li compositions. We will use compositions from 0.1 to 0.5 Li. For the solid structure, we will first create an Al FCC structure, and replace randomly selected atoms with Li. Let’s see how we do this.
structure = pr.create.structure.ase.bulk('Al', cubic=True, a=4.135).repeat(4)
structure.plot3d()
Now we assume we need to create 0.1 composition of Li. Therefore, the number of Li atoms needed are:
comp = 0.1
n_Li = int(comp*len(structure))
n_Li
25
Now we randomly replace 50 Al atoms with Li.
structure[np.random.permutation(len(structure))[:n_Li]] = 'Li'
structure.plot3d()
We can see that some Al atoms are now replaced with Li. We also need to create B32 structures of varying compositions. For that we start with the LiAl B32 structure, and replace randomly selected Li atoms with Al, therby reducing the amount of Li in the structure.
structure = pr.create.structure.ase.read('AlLi_poscar', format='vasp')
structure.plot3d()
Once again, find the number of Li atoms that need to replaced.
n_Li = int((0.5-comp)*len(structure))
n_Li
172
Now replace random Li atoms with Al
rinds = len(structure)//2 + np.random.choice(range(len(structure)//2), n_Li, replace=False)
structure[rinds] = 'Al'
structure.plot3d()
Now we have all the necessary components in place. We can simply create a loop over the compositions and run the calculations. We have prepared the functions which create structures that can be used directly.
for count, comp in tqdm(enumerate([0.1, 0.2, 0.3, 0.4, 0.5])):
structure_fcc, structure_b32 = create_structures(pr, comp, repetitions=4)
job_sol = pr.create.job.Calphy("x%d_sol"%count, delete_aborted_job=True)
job_sol.structure = structure_fcc
job_sol.potential = "LiAl-atomicrex"
job_sol.server.cores = 2
job_sol.calc_free_energy(temperature=[800, 1000],
pressure=0,
reference_phase="solid",
n_equilibration_steps=5000,
n_switching_steps=5000)
job_sol.run()
job_lqd = pr.create.job.Calphy("x%d_lqd"%count, delete_aborted_job=True)
job_lqd.structure = structure_fcc
job_lqd.potential = "LiAl-atomicrex"
job_lqd.server.cores = 2
job_lqd.calc_free_energy(temperature=[800, 1000],
pressure=0,
reference_phase="liquid",
n_equilibration_steps=5000,
n_switching_steps=5000)
job_lqd.run()
job_alli = pr.create.job.Calphy("x%d_alli"%count, delete_aborted_job=True)
job_alli.structure = structure_b32
job_alli.potential = "LiAl-atomicrex"
job_alli.server.cores = 2
job_alli.calc_free_energy(temperature=[800, 1000],
pressure=0,
reference_phase="solid",
n_equilibration_steps=5000,
n_switching_steps=5000)
job_alli.run()
2022-06-07 20:55:03,281 - pyiron_log - WARNING - The job x0_sol is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:03,786 - pyiron_log - WARNING - The job x0_lqd is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:04,280 - pyiron_log - WARNING - The job x0_alli is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:04,873 - pyiron_log - WARNING - The job x1_sol is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:05,356 - pyiron_log - WARNING - The job x1_lqd is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:05,652 - pyiron_log - WARNING - The job x1_alli is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:05,996 - pyiron_log - WARNING - The job x2_sol is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:06,287 - pyiron_log - WARNING - The job x2_lqd is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:06,577 - pyiron_log - WARNING - The job x2_alli is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:06,929 - pyiron_log - WARNING - The job x3_sol is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:07,219 - pyiron_log - WARNING - The job x3_lqd is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:07,506 - pyiron_log - WARNING - The job x3_alli is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:07,815 - pyiron_log - WARNING - The job x4_sol is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:08,107 - pyiron_log - WARNING - The job x4_lqd is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2022-06-07 20:55:08,399 - pyiron_log - WARNING - The job x4_alli is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
Pack the project#
pr.pack("LiAl")
pr = Project("LiAl_analysis")
pr.unpack("LiAl")
pr
{'groups': ['phase_diagram'], 'nodes': []}
Analysing the results#
Now we can analyse the results of the above calculations. First we create an array of the composition values.
comp = np.arange(0, 0.6, 0.1)
For the initial set of analysis, we will choose a temperature of 800 K.
temp = 800
The calculations we ran already has the free energy at all temperatures from 800-1000 K. We need to extract the free energy at the correct temperature. The helpers.py file in the folder contains some helper functions for this notebook. We provide a fe_at method which can extract the free energy at the required temperature. Let’s take a look at the method.
fe_at?
Signature: fe_at(p, temp, threshold=0.1)
Docstring:
Get the free energy at a given temperature
Parameters
----------
p: pyiron Job
Pyiron job with calculated free energy and temperature
temp: float
Required temperature
threshold: optional, default 1E-1
Minimum difference needed between required temperature and temperature found in pyiron job
Returns
-------
float: free energy value at required temperature
File: ~/potentials/phase_diagram/helpers.py
Type: function
For for pure Al calculations, and for each composition, we extract the free energy at 800 K of the FCC, liquid and B32 phases.
fcc = []
b32 = []
lqd = []
fcc.append(fe_at(pr["phase_diagram/xp_sol"], temp))
lqd.append(fe_at(pr["phase_diagram/xp_lqd"], temp))
b32.append(fe_at(pr["phase_diagram/xp_alli"], temp))
for i in range(5):
fcc.append(fe_at(pr["phase_diagram/x%d_sol"%i], temp))
lqd.append(fe_at(pr["phase_diagram/x%d_lqd"%i], temp))
b32.append(fe_at(pr["phase_diagram/x%d_alli"%i], temp))
Plot the results
plt.plot(comp, fcc, '-', color="#e58080")
plt.plot(comp, lqd, '-', color="#66cfff")
plt.plot(comp, b32, '-', color="#ffc766")
plt.plot(comp, fcc, 'o', label='fcc', color="#e58080", markeredgecolor="#424242")
plt.plot(comp, lqd, 'o', label='lqd', color="#66cfff", markeredgecolor="#424242")
plt.plot(comp, b32, 'o', label='b32', color="#ffc766", markeredgecolor="#424242")
plt.xlabel(r"$x_{Li}$")
plt.ylabel(r"F (eV/atom)")
plt.legend()
<matplotlib.legend.Legend at 0x7f6241f3ef10>
Configurational entropy #
In the above example, we had off stoichiometric compositions, but we did not include configurational entropy for the solid structures. The easiest way to do this, which we will use here, is to employ the ideal mixing assumption. In the case of ideal mixing, the configuration entropy of mixing is given by,
We can add this directly to the free energy.
For the liquid phase, addition of configurational entropy explicitely is not needed as is it included in the free energy calculations.
smix = -kb*(comp*np.log(comp) + (1-comp)*np.log(1-comp))
smix[0] = 0
/tmp/ipykernel_1698/1902816561.py:1: RuntimeWarning: divide by zero encountered in log
smix = -kb*(comp*np.log(comp) + (1-comp)*np.log(1-comp))
/tmp/ipykernel_1698/1902816561.py:1: RuntimeWarning: invalid value encountered in multiply
smix = -kb*(comp*np.log(comp) + (1-comp)*np.log(1-comp))
fcc_mix = np.array(fcc)-temp*smix
b32_mix = np.array(b32)-temp*smix
lqd_mix = np.array(lqd)
plt.plot(comp, fcc_mix, '-', color="#e58080")
plt.plot(comp, lqd_mix, '-', color="#66cfff")
plt.plot(comp, b32_mix, '-', color="#ffc766")
plt.plot(comp, fcc_mix, 'o', label='fcc', color="#e58080", markeredgecolor="#424242")
plt.plot(comp, lqd_mix, 'o', label='lqd', color="#66cfff", markeredgecolor="#424242")
plt.plot(comp, b32_mix, 'o', label='b32', color="#ffc766", markeredgecolor="#424242")
plt.xlabel(r"$x_{Li}$")
plt.ylabel(r"F (eV/atom)")
plt.legend()
<matplotlib.legend.Legend at 0x7f62413ba7c0>
To obtain the results in a computationally efficient way, we used composition along every 0.1 Li. We can fit a 3rd order polynomial to these points to get a finer grid. We will use numpy.polyfit for this purpose.
Let’s first define a finer composition grid
comp_grid = np.linspace(0, 0.5, 1000)
Now we fit the free energy values and use this fit to revaluate the free energy on the finer grid.
fcc_fit = np.polyfit(comp, fcc_mix, 3)
fcc_fe = np.polyval(fcc_fit, comp_grid)
lqd_fit = np.polyfit(comp, lqd_mix, 3)
lqd_fe = np.polyval(lqd_fit, comp_grid)
b32_fit = np.polyfit(comp, b32_mix, 3)
b32_fe = np.polyval(b32_fit, comp_grid)
Plot the fits
plt.plot(comp_grid, fcc_fe, '-', color="#e58080")
plt.plot(comp_grid, lqd_fe, '-', color="#66cfff")
plt.plot(comp_grid, b32_fe, '-', color="#ffc766")
plt.plot(comp, fcc_mix, 'o', label='fcc', color="#e58080", markeredgecolor="#424242")
plt.plot(comp, lqd_mix, 'o', label='lqd', color="#66cfff", markeredgecolor="#424242")
plt.plot(comp, b32_mix, 'o', label='b32', color="#ffc766", markeredgecolor="#424242")
plt.xlabel(r"$x_{Li}$")
plt.ylabel(r"F (eV/atom)")
plt.legend()
<matplotlib.legend.Legend at 0x7f6259f01910>
Now in order to identify regions we need to construct common tangents to free energy curves. To make this easier, we convert to free energy of mixing relative to the liquid. For this purpose, we subtract a straight line connecting the free energy of the pure Al liquid to the 0.5 Li liquid, from the calculated free energy. This same line will then be subtracted from the fcc and b32 curves. Our helpers.py module previously introduced includes helper methods to do this.
lqd_fe_norm, slope, intercept = normalise_fe(lqd_fe, comp_grid)
fcc_fe_norm = fcc_fe-(slope*comp_grid + intercept)
b32_fe_norm = b32_fe-(slope*comp_grid + intercept)
plt.plot(comp_grid, fcc_fe_norm, '-', color="#e58080", label='fcc')
plt.plot(comp_grid, lqd_fe_norm, '-', color="#66cfff", label='lqd')
plt.plot(comp_grid, b32_fe_norm, '-', color="#ffc766", label='b32')
plt.xlabel(r"$x_{Li}$")
plt.ylabel(r"F (eV/atom)")
plt.legend()
<matplotlib.legend.Legend at 0x7f622006dc70>
Common tangent constructions #
With the above curves, we can move on to common tangent constructions to identify regions where the two phases coexist. In order to calculate the common tangent for two curves \(f\) and \(g\), we can solve the following set of equationns:
\(x_1\) and \(x_2\) are the endpoints of the common tangent.
The fitting is done using the fsolve method from scipy. Once again, helpers.py module offers a function to do this fitting.
find_common_tangent?
Signature: find_common_tangent(fe1, fe2, guess_range)
Docstring:
Do a common tangent construction between two free energy curves.
Parameters
----------
fe1: numpy array
first free energy curve
fe2: numpy array
second free energy curve
guess_range: list of floats length 2
The guess range to find end points of the common tangent
Returns
-------
res: list of floats length 2
The end points of the common tangent
File: ~/potentials/phase_diagram/helpers.py
Type: function
We will use this method to find the common tangent between the fcc and b32 curves
ct = find_common_tangent(fcc_fit, b32_fit, [0.0, 0.25])
ct
array([0.24357731, 0.32399064])
Note that we provided the polynomials that describe the free energy curves, obtained from fitting to the function. We have obtained \(x_1\) and \(x_2\). We can plot the common tangent now.
plt.plot(comp_grid, fcc_fe_norm, '-', color="#e58080", label='fcc')
plt.plot(comp_grid, lqd_fe_norm, '-', color="#66cfff", label='lqd')
plt.plot(comp_grid, b32_fe_norm, '-', color="#ffc766", label='b32')
plt.plot(ct, [np.polyval(fcc_fit, ct[0])-(slope*ct[0] + intercept),
np.polyval(b32_fit, ct[1])-(slope*ct[1] + intercept)], color="#424242")
plt.xlabel(r"$x_{Li}$")
plt.ylabel(r"F (eV/atom)")
plt.legend()
<matplotlib.legend.Legend at 0x7f62414b4b80>
Lets take a moment to analyse this plot. We obtained 0.24 and 0.32 as the end points of the common tangent. From the figure, this means that FCC is the most stable structure upto \(x_{Li} = 0.24\). Between \(0.24 \leq x_{Li} \leq 0.32\), both fcc and b32 phases coexist. Finally, for \(x_{Li} > 0.32\), b32 is the most stable phase.
Therefore we have obtained already a slice of the phase diagram at 800 K. We can now put all of this methods together and easily calculate for different temperatures.
temp = 900
fcc = []
b32 = []
lqd = []
fcc.append(fe_at(pr["phase_diagram/xp_sol"], temp))
lqd.append(fe_at(pr["phase_diagram/xp_lqd"], temp))
b32.append(fe_at(pr["phase_diagram/xp_alli"], temp))
for i in range(5):
fcc.append(fe_at(pr["phase_diagram/x%d_sol"%i], temp))
lqd.append(fe_at(pr["phase_diagram/x%d_lqd"%i], temp))
b32.append(fe_at(pr["phase_diagram/x%d_alli"%i], temp))
fcc_mix = np.array(fcc)-temp*smix
b32_mix = np.array(b32)-temp*smix
lqd_mix = np.array(lqd)
fcc_fit = np.polyfit(comp, fcc_mix, 3)
fcc_fe = np.polyval(fcc_fit, comp_grid)
lqd_fit = np.polyfit(comp, lqd_mix, 3)
lqd_fe = np.polyval(lqd_fit, comp_grid)
b32_fit = np.polyfit(comp, b32_mix, 3)
b32_fe = np.polyval(b32_fit, comp_grid)
lqd_fe_norm, slope, intercept = normalise_fe(lqd_fe, comp_grid)
fcc_fe_norm = fcc_fe-(slope*comp_grid + intercept)
b32_fe_norm = b32_fe-(slope*comp_grid + intercept)
ct = find_common_tangent(fcc_fit, b32_fit, [0.0, 0.25])
print(ct)
plt.plot(comp_grid, fcc_fe_norm, '-', color="#e58080", label='fcc')
plt.plot(comp_grid, lqd_fe_norm, '-', color="#66cfff", label='lqd')
plt.plot(comp_grid, b32_fe_norm, '-', color="#ffc766", label='b32')
plt.plot(ct, [np.polyval(fcc_fit, ct[0])-(slope*ct[0] + intercept),
np.polyval(b32_fit, ct[1])-(slope*ct[1] + intercept)], color="#424242")
plt.xlabel(r"$x_{Li}$")
plt.ylabel(r"F (eV/atom)")
plt.legend()
[0.25254351 0.3283028 ]
<matplotlib.legend.Legend at 0x7f62413f2dc0>
Now we continue this process for temperatures of 950 K.
temp = 950
fcc = []
b32 = []
lqd = []
fcc.append(fe_at(pr["phase_diagram/xp_sol"], temp))
lqd.append(fe_at(pr["phase_diagram/xp_lqd"], temp))
b32.append(fe_at(pr["phase_diagram/xp_alli"], temp))
for i in range(5):
fcc.append(fe_at(pr["phase_diagram/x%d_sol"%i], temp))
lqd.append(fe_at(pr["phase_diagram/x%d_lqd"%i], temp))
b32.append(fe_at(pr["phase_diagram/x%d_alli"%i], temp))
fcc_mix = np.array(fcc)-temp*smix
b32_mix = np.array(b32)-temp*smix
lqd_mix = np.array(lqd)
fcc_fit = np.polyfit(comp, fcc_mix, 3)
fcc_fe = np.polyval(fcc_fit, comp_grid)
lqd_fit = np.polyfit(comp, lqd_mix, 3)
lqd_fe = np.polyval(lqd_fit, comp_grid)
b32_fit = np.polyfit(comp, b32_mix, 3)
b32_fe = np.polyval(b32_fit, comp_grid)
lqd_fe_norm, slope, intercept = normalise_fe(lqd_fe, comp_grid)
fcc_fe_norm = fcc_fe-(slope*comp_grid + intercept)
b32_fe_norm = b32_fe-(slope*comp_grid + intercept)
plt.plot(comp_grid, fcc_fe_norm, '-', color="#e58080", label='fcc')
plt.plot(comp_grid, lqd_fe_norm, '-', color="#66cfff", label='lqd')
plt.plot(comp_grid, b32_fe_norm, '-', color="#ffc766", label='b32')
ct = find_common_tangent(fcc_fit, lqd_fit, [0.1, 0.3])
print(ct)
ct1 = find_common_tangent(b32_fit, lqd_fit, [0.3, 0.5])
print(ct1)
plt.plot(comp_grid, fcc_fe_norm, '-', color="#e58080", label='fcc')
plt.plot(comp_grid, lqd_fe_norm, '-', color="#66cfff", label='lqd')
plt.plot(comp_grid, b32_fe_norm, '-', color="#ffc766", label='b32')
plt.plot(ct, [np.polyval(fcc_fit, ct[0])-(slope*ct[0] + intercept),
np.polyval(lqd_fit, ct[1])-(slope*ct[1] + intercept)], color="#424242")
plt.plot(ct1, [np.polyval(b32_fit, ct1[0])-(slope*ct1[0] + intercept),
np.polyval(lqd_fit, ct1[1])-(slope*ct1[1] + intercept)], color="#424242")
plt.xlabel(r"$x_{Li}$")
plt.ylabel(r"F (eV/atom)")
plt.legend()
[0.23139167 0.26587513]
[0.35619186 0.31664302]
<matplotlib.legend.Legend at 0x7f62096d3a90>
Now let’s put together all the common tangent constructions to arrive at the phase diagram.
Since we used a very small system, and switching times, the obtained phase diagram will not be accurate. The results shown here uses a system size of about 4000 atoms for each structure, and 50 ps of equilibration and switching time.
Let’s first look at the phase diagram with the calculated points.
And the phase diagram
