Skip to content

Commit

Permalink
Update chemical_potentials code, tests and tutorial
Browse files Browse the repository at this point in the history
  • Loading branch information
@kavanase
kavanase committed Aug 23, 2023
1 parent e5b2e86 commit 4173840
Show file tree
Hide file tree
Showing 3 changed files with 239 additions and 16 deletions.
111 changes: 111 additions & 0 deletions doped/chemical_potentials.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1509,6 +1509,117 @@ def cplap_input(self, dependent_variable=None, filename="input.dat"):
print(int(v), k, end=" ")
print(f"{i['formation_energy']} # num_atoms, element, formation_energy")

def to_LaTeX_table(self, columns=1, sort_by_energy=False, prune_polymorphs=True):
"""
A very simple function to print out the competing phase formation
energies in a LaTeX table format. Needs the mhchem package to work and
does *not* use the booktabs package -- change hline to toprule, midrule
and bottomrule if you want to use booktabs style.
Args:
columns (int):
Number of columns to print out the table in; either 1 or 2
sort_by_energy (bool):
If True, sorts the table by formation energy (highest to lowest).
Default sorting by formula.
prune_polymorphs (bool):
Whether to only print out the lowest energy polymorphs for each composition.
Default is True.
Returns:
str: LaTeX table string
"""
if columns not in [1, 2]:
raise ValueError("columns must be either 1 or 2")
# done in the pyscfermi report style
data = copy.deepcopy(self.data)

if any("kpoints" not in item for item in data):
raise (
ValueError(
"kpoints need to be present in data; run CompetingPhasesAnalyzer.from_vaspruns "
"instead of from_csv"
)
)

if prune_polymorphs: # only keep the lowest energy polymorphs
formation_energy_df = _calculate_formation_energies(data, self.elemental_energies)
indices = formation_energy_df.groupby("formula")["energy_per_atom"].idxmin()
pruned_df = formation_energy_df.loc[indices]
data = pruned_df.to_dict(orient="records")

if sort_by_energy:
data = sorted(data, key=lambda x: x["formation_energy"], reverse=True)
# moves the bulk composition to the top of the list
_move_dict_to_start(data, "formula", self.bulk_composition.reduced_formula)

string = "\\begin{table}[h]\n\\centering\n"
string += (
"\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{"
+ self.bulk_composition.reduced_formula
+ "} and all competing phases, with k-meshes used in calculations."
+ ("}\n" if not prune_polymorphs else " Only the lowest energy polymorphs are included}\n")
)
string += "\\label{tab:competing_phase_formation_energies}\n"
if columns == 1:
string += "\\begin{tabular}{ccc}\n"
string += "\\hline\n"
string += "Formula & k-mesh & $\\Delta E_f$ (eV) \\\\ \\hline \n"
for i in data:
kpoints = i["kpoints"].split("x")
fe = i["formation_energy"]
string += (
"\\ce{"
+ i["formula"]
+ "} & "
+ f"{kpoints[0]}$\\times${kpoints[1]}$\\times${kpoints[2]}"
" & " + f"{fe:.3f} \\\\ \n"
)

elif columns == 2:
string += "\\begin{tabular}{ccc|ccc}\n"
string += "\\hline\n"
string += (
"Formula & k-mesh & $\\Delta E_f$ (eV) & Formula & k-mesh & $\\Delta E_f$ (eV)\\\\ "
"\\hline \n"
)

mid = len(data) // 2
first_half = data[:mid]
last_half = data[mid:]

for i, j in zip(first_half, last_half):
kpoints = i["kpoints"].split("x")
fe = i["formation_energy"]
kpoints2 = j["kpoints"].split("x")
fe2 = j["formation_energy"]
string += (
"\\ce{"
+ i["formula"]
+ "} & "
+ f"{kpoints[0]}$\\times${kpoints[1]}$\\times${kpoints[2]}"
" & "
+ f"{fe:.3f} & "
+ "\\ce{"
+ j["formula"]
+ "} & "
+ f"{kpoints2[0]}$\\times${kpoints2[1]}$\\times${kpoints2[2]}"
" & " + f"{fe2:.3f} \\\\ \n"
)

string += "\\hline\n"
string += "\\end{tabular}\n"
string += "\\end{table}"

return string


def _move_dict_to_start(data, key, value):
for index, item in enumerate(data):
if key in item and item[key] == value:
data.insert(0, data.pop(index))
return


def combine_extrinsic(first, second, extrinsic_species):
# TODO: Can we just integrate this to `CompetingPhaseAnalyzer`, so you just pass in a list of
Expand Down
103 changes: 88 additions & 15 deletions examples/dope_chemical_potentials.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -377,7 +377,9 @@
"\n",
"You may want to change the default `ENCUT` (350 eV) or k-point densities that the convergence tests span (5 - 60 kpoints/Å<sup>3</sup> for semiconductors & insulators and 40 - 120 kpoints/Å<sup>3</sup> for metals in steps of 5 kpoints/Å<sup>3</sup>). Note that `ISMEAR = -5` is used for metals by default (better kpoint convergence for energies but should not be used during metal geometry relaxations) and k-point convergence testing is not required for molecules (Γ-point sampling is sufficient).\n",
"\n",
"Note that `doped` generates \"molecule in a box\" structures for the gaseous elemental phases H<sub>2</sub>, O<sub>2</sub>, N<sub>2</sub>, F<sub>2</sub> and Cl<sub>2</sub>. The molecule is placed in a 30 Å x 30 Å x 30 Å box, and relaxed with Γ-point-only k-point sampling.\n",
"Note that `doped` generates \"molecule in a box\" structures for the gaseous elemental phases\n",
"H<sub>2</sub>, O<sub>2</sub>, N<sub>2</sub>, F<sub>2</sub> and Cl<sub>2</sub>. The molecule is placed in\n",
" a slightly-symmetry-broken (to avoid metastable electronic solutions) 30 Å cuboid box, and relaxed with Γ-point-only k-point sampling.\n",
"\n",
"The kpoints convergence calculations are set up with:"
]
Expand Down Expand Up @@ -525,7 +527,11 @@
"id": "7e57cb66",
"metadata": {},
"source": [
"Remember that the final `ENCUT` used for the energy calculations should be the same as for your host material & defects, and that you may still need to account for Pulay stress by increasing `ENCUT` for the geometry relaxations (usual rule of thumb being 1.3*converged `ENCUT`) or re-relaxing each structure until the volume change is minimal (roughly <0.3%)."
"Remember that the final `ENCUT` used for the energy calculations should be the same as for your host\n",
"material & defects, and that you may still need to account for Pulay stress by increasing `ENCUT` for\n",
"the geometry relaxations (a typical rule of thumb being 1.3*converged `ENCUT`) or re-relaxing each\n",
"structure until the volume change is minimal (roughly <0.3%). This is not the case for the\n",
"molecule-in-a-box competing phases however, due to the large simulation box size and fixed volume."
]
},
{
Expand Down Expand Up @@ -730,16 +736,16 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 1,
"id": "e273f8fa",
"metadata": {
"ExecuteTime": {
"end_time": "2023-08-21T11:54:14.575272Z",
"start_time": "2023-08-21T11:54:14.548811Z"
},
"collapsed": false,
"jupyter": {
"outputs_hidden": false
},
"ExecuteTime": {
"end_time": "2023-08-23T21:18:44.960159Z",
"start_time": "2023-08-23T21:18:39.250023Z"
}
},
"outputs": [],
Expand All @@ -750,12 +756,12 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 2,
"id": "dd65cb61-bd10-41c8-bf6d-15398841f226",
"metadata": {
"ExecuteTime": {
"end_time": "2023-08-21T11:54:16.335532Z",
"start_time": "2023-08-21T11:54:15.916656Z"
"end_time": "2023-08-23T21:18:45.128309Z",
"start_time": "2023-08-23T21:18:44.961171Z"
}
},
"outputs": [
Expand Down Expand Up @@ -789,16 +795,16 @@
},
{
"cell_type": "code",
"execution_count": 23,
"execution_count": 3,
"id": "a3d1b363",
"metadata": {
"ExecuteTime": {
"end_time": "2023-08-21T11:51:09.266661Z",
"start_time": "2023-08-21T11:51:09.052205Z"
},
"collapsed": false,
"jupyter": {
"outputs_hidden": false
},
"ExecuteTime": {
"end_time": "2023-08-23T21:18:45.377891Z",
"start_time": "2023-08-23T21:18:45.129641Z"
}
},
"outputs": [
Expand Down Expand Up @@ -844,6 +850,73 @@
"cpa.from_vaspruns(all_paths)"
]
},
{
"cell_type": "markdown",
"source": [
"### Print Formation Energies as a `LaTeX` Table\n",
"You may want to print out the parsed calculation results for your competing phases as a LaTeX table, to\n",
"be included in supplementary information of journal articles or in theses, to aid open-science and reproducibility.\n",
"The formation energies are automatically saved to a `csv` file, which can be converted to LaTeX (see\n",
"https://www.tablesgenerator.com/latex_tables) or\n",
"Word formats, but you can also print out a nicely-formatted LaTeX table using the `to_LaTeX_table()`\n",
"method as shown belo –\n",
"either as a one- or two-column table; ordered alphabetially or in descending order of energies."
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
"execution_count": 4,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\\begin{table}[h]\n",
"\\centering\n",
"\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{ZrO2} and all competing phases, with k-meshes used in calculations. Only the lowest energy polymorphs are included}\n",
"\\label{tab:competing_phase_formation_energies}\n",
"\\begin{tabular}{ccc}\n",
"\\hline\n",
"Formula & k-mesh & $\\Delta E_f$ (eV) \\\\ \\hline \n",
"\\ce{ZrO2} & 3$\\times$3$\\times$3 & -10.975 \\\\ \n",
"\\ce{O2} & 2$\\times$2$\\times$2 & 0.000 \\\\ \n",
"\\ce{Zr} & 9$\\times$9$\\times$5 & 0.000 \\\\ \n",
"\\ce{Zr2O} & 5$\\times$5$\\times$2 & -5.729 \\\\ \n",
"\\ce{Zr3O} & 5$\\times$5$\\times$5 & -5.987 \\\\ \n",
"\\hline\n",
"\\end{tabular}\n",
"\\end{table}\n"
]
}
],
"source": [
"table = cpa.to_LaTeX_table()\n",
"print(table)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-08-23T21:18:45.385574Z",
"start_time": "2023-08-23T21:18:45.382991Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"```{note}\n",
"The table does not use `booktabs` by default but you can change the horizontal lines (`hline`) to\n",
"toprule, midrule and bottomrule to get the desired style effect. The output also relies on using the\n",
"`\\ce{}` command from the `mhchem` LaTeX package to get the correct formatting for the chemical formulae.\n",
"```"
],
"metadata": {
"collapsed": false
}
},
{
"attachments": {},
"cell_type": "markdown",
Expand Down
41 changes: 40 additions & 1 deletion tests/test_chemical_potentials.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,7 +217,46 @@ def test_cplap_input(self):
assert (
"1 Zr 2 O -10.951109995000005\n" not in contents
) # shouldn't be included as is a higher energy polymorph of the bulk phase
assert not any("2 Zr 1 O" in i for i in contents) # non-bordering phase
assert all("2 Zr 1 O" not in i for i in contents) # non-bordering phase

def test_latex_table(self):
cpa = chemical_potentials.CompetingPhasesAnalyzer(self.stable_system)
path = self.path / "ZrO2"
cpa.from_vaspruns(path=path, folder="relax", csv_fname=self.csv_path)

with self.assertRaises(ValueError):
cpa.to_LaTeX_table(columns="M")

string = cpa.to_LaTeX_table(columns=1)
assert (
string[28:209]
== "\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{ZrO2} and all "
"competing phases, with k-meshes used in calculations. Only the lowest energy polymorphs "
"are included"
)
assert len(string) == 586
assert string.split("hline")[1] == "\nFormula & k-mesh & $\\Delta E_f$ (eV) \\\\ \\"
assert string.split("hline")[2][2:45] == "\\ce{ZrO2} & 3$\\times$3$\\times$3 & -10.975 \\"

string = cpa.to_LaTeX_table(columns=2)
assert (
string[28:209]
== "\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{ZrO2} and all "
"competing phases, with k-meshes used in calculations. Only the lowest energy polymorphs "
"are included"
)
assert (
string.split("hline")[1]
== "\nFormula & k-mesh & $\\Delta E_f$ (eV) & Formula & k-mesh & $\\Delta E_f$ (eV)\\\\ \\"
)

assert string.split("hline")[2][2:45] == "\\ce{ZrO2} & 3$\\times$3$\\times$3 & -10.975 &"
assert len(string) == 579

cpa_csv = chemical_potentials.CompetingPhasesAnalyzer(self.stable_system)
cpa_csv.from_csv(self.csv_path)
with self.assertRaises(ValueError):
cpa_csv.to_LaTeX_table(columns=2)


class BoxedMoleculesTestCase(unittest.TestCase):
Expand Down

0 comments on commit 4173840

Please sign in to comment.