Running Materials Project Compatible workflows using VASP / DFT¶
As of early 2025, the majority of the calculations in the Materials Project (MP) database are performed with PBE GGA and PBE+U density funcitonal theory (DFT) calculations, with about a quarter coming from r2SCAN calculations instead.
In this example, you will use the Materials Project workflows and corrections/mixing schemes for PBE and PBE+U, and for r2SCAN, to generate a partial compositional phase diagram for the V-O compositional space. The mixing and corrections schemes are used to better align PBE, PBE+U, and r2SCAN convex hulls, and to correct known errors in approximate DFT to predict more accurate enthalpies of formation.
Note: The Materials Project is currently recomputing the entirety of its database with updated r2SCAN workflows. We will focus first on the current set of calculations, and mention these updated workflows at the end of the notebook.
Let’s start by investigating the input sets used by MP
from __future__ import annotations
from pymatgen.io.vasp.sets import MPRelaxSet, MPScanRelaxSet
print( # noqa: T201
"PBE GGA / +U INCAR:\n"
+ "\n".join(f" {k}: {v}" for k, v in MPRelaxSet.CONFIG["INCAR"].items())
)
print( # noqa: T201
"r2SCAN meta-GGA INCAR:\n"
+ "\n".join(f" {k}: {v}" for k, v in MPScanRelaxSet.CONFIG["INCAR"].items())
)
While you can set up MP compatible VASP calculations by hand, it’s usually best practice to use the pre-defined workflows in atomate2.
For PBE GGA/+U calculations, the workflow is atomate2.vasp.flows.mp.MPGGADoubleRelaxStaticMaker. For r2SCAN, the workflow is atomate2.vasp.flows.mp.MPMetaGGADoubleRelaxStaticMaker.
We start by pulling a few structures from the Materials Project API in the V-O system, importantly including the unaries.
from mp_api.client import MPRester
with MPRester() as mpr:
docs = mpr.materials.summary.search(
material_ids=["mp-146", "mp-12957", "mp-18937", "mp-541404"],
fields=["material_id", "structure"],
)
structures = {doc.material_id.string: doc.structure.to_primitive() for doc in docs}
for mp_id, structure in structures.items():
print(structure.to(fmt="poscar", comment=f"{mp_id} {structure.formula}")) # noqa: T201
If you don’t have access to VASP, or just want to run the tutorial example:¶
All test data is provided with this example (including raw VASP input and output), just run the following code block:
from mock_vasp import TEST_DIR, mock_vasp
from monty.serialization import loadfn
tasks = loadfn(
TEST_DIR
/ "vasp"
/ "mp_phase_diagram_tutorial"
/ "MP-example-V-O-phase-diagram.json.gz"
)
structures = {
mp_id: task_block["PBEsol Structure Optimization"].input.structure
for mp_id, task_block in tasks.items()
}
Now let’s actually set up the workflows! For these calculations, we will adjust “parallelization settings” to ensure that our calculations run quickly. If you are not familiar with these settings in VASP, it’s best to read the manual section in VASP.
If you run VASP on CPU nodes on a cluster/HPC, you want to modify:
NCORE: how many individual cores work on a single orbital
KPAR: how many k points are treated in parallel
If you run VASP on GPU nodes on a cluster/HPC, you will modify:
NSIM: parallelization over bands for certain optimizers
KPAR
NCORE = 1
To do the latter part, we will use the powerups in atomate2. We will also save our results to JSON-format output using monty.
For parallelization settings, I’m running these on GPU nodes, so I set NCORE = 1, KPAR = 8 (using 8 individual GPUs), and setting NSIM=8 as a reasonable value that does not increase the memory demands too much.’
Note that we also update the names of the jobs and their metadata to make it easier to query output later when we build the phase diagram.
from jobflow import Flow, JobStore, run_locally
from maggma.stores import MemoryStore
from mock_vasp import TEST_DIR
from monty.serialization import dumpfn
from atomate2.vasp.flows.mp import (
MPGGADoubleRelaxStaticMaker,
MPMetaGGADoubleRelaxStaticMaker,
)
from atomate2.vasp.powerups import update_user_incar_settings
jobs = []
for mp_id, structure in structures.items():
pbe_job = MPGGADoubleRelaxStaticMaker().make(structure)
pbe_job.append_name(f" {mp_id}")
pbe_job.update_metadata({"mp_id": mp_id})
r2scan_job = MPMetaGGADoubleRelaxStaticMaker().make(structure)
r2scan_job.append_name(f" {mp_id}")
pbe_job.update_metadata({"mp_id": mp_id})
jobs.extend([pbe_job, r2scan_job])
all_jobs = Flow(jobs)
all_jobs = update_user_incar_settings(all_jobs, {"NCORE": 1, "KPAR": 8, "NSIM": 8})
test_data_dir = TEST_DIR / "vasp" / "mp_phase_diagram_tutorial"
ref_paths = {
k: test_data_dir / v
for k, v in {
"MP GGA static mp-12957": "MP_GGA_static_mp-12957.json.lzma",
"MP GGA relax 2 mp-18937": "MP_GGA_relax_2_mp-18937.json.lzma",
"MP GGA relax 1 mp-18937": "MP_GGA_relax_1_mp-18937.json.lzma",
"MP pre-relax 1 mp-146": "MP_pre-relax_1_mp-146.json.lzma",
"MP meta-GGA relax 2 mp-146": "MP_meta-GGA_relax_2_mp-146.json.lzma",
"MP pre-relax 1 mp-541404": "MP_pre-relax_1_mp-541404.json.lzma",
"MP meta-GGA static mp-146": "MP_meta-GGA_static_mp-146.json.lzma",
"MP GGA static mp-146": "MP_GGA_static_mp-146.json.lzma",
"MP meta-GGA static mp-12957": "MP_meta-GGA_static_mp-12957.json.lzma",
"MP meta-GGA static mp-541404": "MP_meta-GGA_static_mp-541404.json.lzma",
"MP meta-GGA relax 2 mp-541404": "MP_meta-GGA_relax_2_mp-541404.json.lzma",
"MP GGA relax 2 mp-146": "MP_GGA_relax_2_mp-146.json.lzma",
"MP meta-GGA relax 2 mp-12957": "MP_meta-GGA_relax_2_mp-12957.json.lzma",
"MP pre-relax 1 mp-12957": "MP_pre-relax_1_mp-12957.json.lzma",
"MP GGA relax 1 mp-12957": "MP_GGA_relax_1_mp-12957.json.lzma",
"MP GGA static mp-541404": "MP_GGA_static_mp-541404.json.lzma",
"MP GGA relax 1 mp-146": "MP_GGA_relax_1_mp-146.json.lzma",
"MP GGA static mp-18937": "MP_GGA_static_mp-18937.json.lzma",
"MP GGA relax 2 mp-12957": "MP_GGA_relax_2_mp-12957.json.lzma",
"MP pre-relax 1 mp-18937": "MP_pre-relax_1_mp-18937.json.lzma",
"MP meta-GGA relax 2 mp-18937": "MP_meta-GGA_relax_2_mp-18937.json.lzma",
"MP GGA relax 1 mp-541404": "MP_GGA_relax_1_mp-541404.json.lzma",
"MP meta-GGA static mp-18937": "MP_meta-GGA_static_mp-18937.json.lzma",
"MP GGA relax 2 mp-541404": "MP_GGA_relax_2_mp-541404.json.lzma",
}.items()
}
job_store = JobStore(MemoryStore(), additional_stores={"data": MemoryStore()})
with mock_vasp(ref_paths) as vasp:
responses = run_locally(
all_jobs,
store=job_store,
create_folders=True,
ensure_success=True,
)
If you ran this locally, you can collect the output of the jobs from the local directory structure. We use the formula of the original structures and calc_type (name of the functional used plus the type of job performed, like a static or a relaxation) to distinguish jobs. We will use tasks to plot phase diagrams.
from pathlib import Path
from emmet.core.tasks import TaskDoc
formula_to_mpid = {s.formula: mpid for mpid, s in structures.items()}
task_list = [TaskDoc.from_directory(p) for p in Path(".").glob("job_*") if p.is_dir()]
tasks = {mpid: {} for mpid in structures}
for task in task_list:
mpid = formula_to_mpid.get(task.structure.formula)
ct = task.calc_type.value
if ct in tasks[mpid]:
nmatch = len([k for k in tasks[mpid] if ct in k])
ct += f" {nmatch + 1}"
tasks[mpid][ct] = task
dumpfn(tasks, "V-O_MP_task_docs.json.gz")
Alternatively, you could use the following code to obtain the output using the job names directly. It’s a good exercise to see how you can use all_output in the following phase diagram plots.
from atomate2.utils.testing import get_job_uuid_name_map
job_uuid_to_name = get_job_uuid_name_map(all_jobs)
all_output = {
name: responses[uuid][1].output for uuid, name in job_uuid_to_name.items()
}
Last, and this is the way you should go in high-throughput, we can query the job store for output. The contents of the job store are “deserialized,” so we reconstruct them with their document model, emmet.core.tasks.TaskDoc.
Note that we only need the static calculations (with task_type of GGA Static, GGA+U Static, and r2SCAN Static) to build the phase diagram, as these most accurately capture the total energy. That’s because of various electronic structure calculation reasons (e.g., as the lattice changes during a geomtry optimization, the k point grid may no longer have high-enough density, or it may represent the wrong lattice symmetry).
from emmet.core.tasks import TaskDoc
all_output = {}
for mp_id in structures:
all_output[mp_id] = {}
for entry in job_store.query(
{"metadata.mp_id": mp_id, "name": {"$regex": "static"}}
):
all_output[mp_id][entry["output"]["calc_type"]] = TaskDoc(**entry["output"])
We will use these static calculations to build a compositional phase diagram. First, we have to apply mixing and correction scheme corrections to these tasks. These are applied by pymatgen.entries.mixing_scheme.MaterialsProjectDFTMixingScheme - we only have to pass a small representation of each calculation called a ComputedStructureEntry, to the mixer, and it applies all empirical corrections / functional mixing adjustments.
Note: MaterialsProjectDFTMixingScheme will give sensible output only if you have complete phase space information (elemental unary compounds and relevant intermediates). We will cover the case of applying corrections without adjusting the PBE and r2SCAN hulls in a moment.
To best demonstrate how the mixing corrections work, we will purposefully omit one of the r2SCAN calculations.
from pymatgen.entries.mixing_scheme import MaterialsProjectDFTMixingScheme
static_entries = []
for mpid, calcs in all_output.items():
for calc_type, task_doc in calcs.items():
if mpid == "mp-18937" and calc_type.lower() == "r2scan static":
continue
if "static" in calc_type.lower():
cse = task_doc.structure_entry
cse.entry_id = f"{cse.parameters['run_type']} {mpid}"
cse.data["material_id"] = mpid
static_entries.append(cse)
mixer = MaterialsProjectDFTMixingScheme()
corrected_entries = mixer.process_entries(static_entries)
for ce in corrected_entries:
print(ce.entry_id, ce.energy_adjustments) # noqa: T201
Now let’s actually plot the phase diagram! Note that the one calculation where we purposefully skipped the r2SCAN calculation, mp-18937, is represented by a PBE+U calculation. You can see this by hovering your cursor over that particular data point.
Try to compare your phase diagram to the one that MP calculates.
What happens if you do not skip this calculation? How does the convex hull change?
from pymatgen.analysis.phase_diagram import PhaseDiagram
pd = PhaseDiagram(corrected_entries)
fig = pd.get_plot()
fig.update_layout(width=1000, height=700, plot_bgcolor="white", paper_bgcolor="white")
To apply just corrections/preliminary mixing to the PBE and PBE+U entries, you would use pymatgen.entries.compatibility.MaterialsProject2020Compatibility. This is useful if you want to compare to older MP data in the legacy site. The correction scheme will raise a CompatibilityError if you attempt to apply it to non-PBE/+U entries.
from pymatgen.entries.compatibility import MaterialsProject2020Compatibility
from pymatgen.entries.computed_entries import ComputedStructureEntry
compat_scheme = MaterialsProject2020Compatibility()
gga_gga_u_corrected_entries = []
for mpid, calcs in tasks.items():
for calc_type, task_doc in calcs.items():
if all(kwd in calc_type.lower() for kwd in ("gga", "static")):
cse_dict = task_doc.structure_entry.as_dict()
cse_dict["energy_adjustments"] = compat_scheme.get_adjustments(
task_doc.structure_entry
)
corrected_cse = ComputedStructureEntry.from_dict(cse_dict)
print( # noqa: T201
mpid, corrected_cse.formula, calc_type, corrected_cse.energy_adjustments
)
Note on revised MP input sets / workflows¶
The Materials Project recently undertook a significant benchmarking effort to verify the accuracy of its DFT calculations. That ultimately led to the construction of new r2SCAN input sets, which are currently being used to recompute all entries in MP.
While these data aren’t currently publicly available, if you want to use these newer input sets / workflows, we list them below:
from pymatgen.io.vasp.sets import MP24RelaxSet
from atomate2.vasp.flows.mp import MP24DoubleRelaxStaticMaker
Then you would use MP24DoubleRelaxStaticMaker in place of MPMetaGGARelaxMaker.