Hello there,
I am trying to reproduce the convex hull represented in Figure 2 of this paper (https://pubs.acs.org/doi/10.1021/acsaem.4c01113).
Right now, I have tried several approaches. I used the supercell program to generate multiple structures at each concentration and also employed the ICET algorithm to produce several structures. I assume that I need to obtain the relaxation energy for each generated configuration, which I am doing using the orbml potential.
My problem arises when I consider the formula for the formation energy, which should be something like:
E_{\text{formation}} = E_{\text{config interest}} - xE_{\text{LiCoO}_2} - (1-x)E_{\text{CoO}_2}
The issue is that this energy should be normalized to formula units. However, for the configurations I generate, I either cannot determine the correct number of formula units or do not fully understand the concept.
My question is: If I have a generated defect configuration, for example, Li_{0.5}Be_{0.5}CoO_2, where I use beryllium as a vacancy placeholder, how can I determine the number of reduced formula units for this configuration, and calculate a convex hull?
Best,
The trick here is to normalize your energy axis by the number of non-Li atoms. Check out grand potential phase diagrams for reference. Here is where the normalization happens.
In grand potential phase diagrams, the axes are normalized grand potential \Phi = G-\mu*N_{Li} vs. x_{Li}, where \Phi in that equation gets divided by number of non-Li atoms before you calculate the convex hull.
For the normalization, it’s easiest to do the energy calculation for the whole cell first and then divide by the number of non-Li atoms. For example, if your structure had 11 atoms and two of them were lithium, you’d take your total free energy (or approximated as just energy) of the 11-atom cell, subtract 2*\mu_{Li}, and then divide by 9 non-Li atoms. Do that for every structure and then calculate the convex hull against Li concentration (i.e., x_{Li} ranging from 0 to 1). If you use pymatgen, it will do this all for you, as seen in the link above.
Thank you very much for your prompt response! I would greatly appreciate it if you could help me with the workflow. Here’s what I have so far:
- Select the stable phases of LiCoO_2 and CoO_2 from the Materials Project.
- Use ICET, Supercell, or another program to generate random configurations for each x in Li_xCoO_2 with x < 1. In this case, I use Be atoms to represent vacancies, following an ICET example.
- Relax each configuration and obtain the free energy.
At this point, I have some doubts:
- Should I remove the Be atoms representing vacancies before relaxation?
- How should I proceed with pymatgen? You pointed me to
GrandPotPDEntry(), so to instantiate this class, I assume I first need to create PDEntry objects for each of my structures. Assuming this is done as follows:
from pymatgen.analysis.phase_diagram import PDEntry
entry = PDEntry(comp, energy, name=structure.composition.reduced_formula)
- Should the
energy parameter be the total energy of the cell or the normalized energy you mentioned?
- Additionally, I know that
structure.composition.reduced_formula in pymatgen will not return a formula like Li_0.5CoO_2. Does this affect the calculations?
Now, when instantiating GrandPotPDEntry(PDEntry):
grand_entry = GrandPotPDEntry(pd_entry, chempot_dict)
- Where do I get the correct chemical potential values?
- What are the free species in this context?
- How can I plot my convex hull from here?
Thank you very much in advance!
Best regards,