Repository: Quantum-Accelerators/quacc (python)

Issue:
Use MP VASP input sets directly without copying them in quacc
### What new feature would you like to see?

This will ensure 1:1 compatibility and reduce the need to maintain them in `presets`. See https://github.com/materialsproject/pymatgen/pull/3484

Pull request description:
Support Pymatgen-style VASP input sets
Closes #1311.

Requires: https://github.com/materialsproject/pymatgen/pull/3615

This PR allows support for Pymatgen-formatted input sets in the `Vasp` calculator and, by extension, all VASP recipes. The `Vasp` calculator takes a new keyword argument `pmg_input_set: DictSet | None = None` that, if specified, would act in the same was a preset by providing a list of default parameters to use. Any class inheriting from [`pymatgen.io.vasp.sets.DictSet`](https://github.com/materialsproject/pymatgen/blob/main/pymatgen/io/vasp/sets.py) is usable. This also means we can use the MP-compatible, pymatgen input sets directly without having to try and duplicate them 1:1 in `quacc`, such as by doing `pmg_input_set = MPRelaxSet`.

I also renamed `mp_relax_job` to `mp_metagga_relax_job` (and similar for the rest) to allow for a future addition of the GGA MP workflows.

Requires: A new release of pymatgen to be minted (bumped in #1622).

Example:

```python
from ase.build import bulk
from quacc.calculators.vasp import Vasp
from pymatgen.io.vasp.sets import MPRelaxSet

atoms = bulk("Fe") * (2, 1, 1)
atoms[0].symbol = "O"
calc = Vasp(atoms, pmg_input_set=MPRelaxSet)
print(calc.parameters)
```

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 99d4645e..3e3d7c9c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -493,7 +493,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 ### Fixed
 
 - VASP: Fixed a scenario where SIGMA was set to 0.05 if the user provided SIGMA > 0.05. This was supposed to happen when ISMEAR = 0 but was previously done when ISMEAR = -5, which is not influenced by SIGMA.
-- Fixed `mp_relax_flow` to use a `#!Python @flow` decorator instead of `#!Python @job`
+- Fixed `mp_metagga_relax_flow` to use a `#!Python @flow` decorator instead of `#!Python @job`
 - Fixed the Q-Chem Quasi-IRC job to run as one `#!Python @job` instead of two.
 
 ## [0.2.3]
diff --git a/docs/dev/recipes/flows.md b/docs/dev/recipes/flows.md
index e90675d0..b3e6d164 100644
--- a/docs/dev/recipes/flows.md
+++ b/docs/dev/recipes/flows.md
@@ -12,7 +12,7 @@ A flow is a Python function decorated by `#!Python @flow` that contains a collec
 
 ### Production Example
 
-A simple, representative flow can be found in [quacc.recipes.vasp.mp.mp_relax_flow][].
+A simple, representative flow can be found in [quacc.recipes.vasp.mp.mp_metagga_relax_flow][].
 
 !!! Note
 
diff --git a/docs/user/recipes/recipes_list.md b/docs/user/recipes/recipes_list.md
index d3292eb7..39f25109 100644
--- a/docs/user/recipes/recipes_list.md
+++ b/docs/user/recipes/recipes_list.md
@@ -229,18 +229,19 @@ The list of available quacc recipes is shown below. The "Req'd Extras" column sp
 
 <center>
 
-| Name                    | Decorator        | Documentation                                   | Req'd Extras |
-| ----------------------- | ---------------- | ----------------------------------------------- | ------------ |
-| VASP Static             | `#!Python @job`  | [quacc.recipes.vasp.core.static_job][]          |              |
-| VASP Relax              | `#!Python @job`  | [quacc.recipes.vasp.core.relax_job][]           |              |
-| VASP Double Relax       | `#!Python @job`  | [quacc.recipes.vasp.core.double_relax_job][]    |              |
-| VASP Slab Static        | `#!Python @job`  | [quacc.recipes.vasp.slabs.static_job][]         |              |
-| VASP Slab Relax         | `#!Python @job`  | [quacc.recipes.vasp.slabs.relax_job][]          |              |
-| VASP Bulk to Slabs      | `#!Python @flow` | [quacc.recipes.vasp.slabs.bulk_to_slabs_flow][] |              |
-| VASP Slab to Adsorbates | `#!Python @flow` | [quacc.recipes.vasp.slabs.slab_to_ads_flow][]   |              |
-| VASP MP Prerelax        | `#!Python @job`  | [quacc.recipes.vasp.mp.mp_relax_job][]          |              |
-| VASP MP Relax           | `#!Python @job`  | [quacc.recipes.vasp.mp.mp_relax_job][]          |              |
-| VASP MP Relax Workflow  | `#!Python @flow` | [quacc.recipes.vasp.mp.mp_relax_flow][]         |              |
-| VASP QMOF Relax         | `#!Python @job`  | [quacc.recipes.vasp.qmof.qmof_relax_job][]      |              |
+| Name                            | Decorator        | Documentation                                   | Req'd Extras |
+| ------------------------------- | ---------------- | ----------------------------------------------- | ------------ |
+| VASP Static                     | `#!Python @job`  | [quacc.recipes.vasp.core.static_job][]          |              |
+| VASP Relax                      | `#!Python @job`  | [quacc.recipes.vasp.core.relax_job][]           |              |
+| VASP Double Relax               | `#!Python @job`  | [quacc.recipes.vasp.core.double_relax_job][]    |              |
+| VASP Slab Static                | `#!Python @job`  | [quacc.recipes.vasp.slabs.static_job][]         |              |
+| VASP Slab Relax                 | `#!Python @job`  | [quacc.recipes.vasp.slabs.relax_job][]          |              |
+| VASP Bulk to Slabs              | `#!Python @flow` | [quacc.recipes.vasp.slabs.bulk_to_slabs_flow][] |              |
+| VASP Slab to Adsorbates         | `#!Python @flow` | [quacc.recipes.vasp.slabs.slab_to_ads_flow][]   |              |
+| VASP MP Meta-GGA Prerelax       | `#!Python @job`  | [quacc.recipes.vasp.mp.mp_metagga_relax_job][]  |              |
+| VASP MP Meta-GGA Relax          | `#!Python @job`  | [quacc.recipes.vasp.mp.mp_metagga_relax_job][]  |              |
+| VASP MP Meta-GGA Static         | `#!Python @job`  | [quacc.recipes.vasp.mp.mp_metagga_static_job][] |              |
+| VASP MP Meta-GGA Relax Workflow | `#!Python @flow` | [quacc.recipes.vasp.mp.mp_metagga_relax_flow][] |              |
+| VASP QMOF Relax                 | `#!Python @job`  | [quacc.recipes.vasp.qmof.qmof_relax_job][]      |              |
 
 </center>
diff --git a/src/quacc/calculators/vasp/params.py b/src/quacc/calculators/vasp/params.py
index 00404ee5..362a120d 100644
--- a/src/quacc/calculators/vasp/params.py
+++ b/src/quacc/calculators/vasp/params.py
@@ -7,6 +7,7 @@ from typing import TYPE_CHECKING
 
 import numpy as np
 from ase.calculators.vasp import Vasp as Vasp_
+from pymatgen.io.ase import AseAtomsAdaptor
 
 from quacc.atoms.core import check_is_metal
 from quacc.utils.kpts import convert_pmg_kpts
@@ -15,6 +16,7 @@ if TYPE_CHECKING:
     from typing import Any, Literal
 
     from ase.atoms import Atoms
+    from pymatgen.io.vasp.sets import DictSet
 
     from quacc.utils.kpts import PmgKpts
 
@@ -349,3 +351,49 @@ def set_pmg_kpts(
         user_calc_params["gamma"] = gamma
 
     return user_calc_params
+
+
+def get_pmg_input_set_params(dict_set: DictSet, atoms: Atoms) -> tuple[dict, Atoms]:
+    """
+    Convert a Pymatgen VASP input set into an ASE-compatible set of
+    calculator parameters.
+
+    Parameters
+    ----------
+    dict_set
+        The Pymatgen VASP input set.
+    atoms
+        The input atoms.
+
+    Returns
+    -------
+    dict
+        The ASE-compatible set of calculator parameters.
+    Atoms
+        The input atoms to match the pymatgen input set.
+    """
+    structure = AseAtomsAdaptor.get_structure(atoms)
+    pmg_input_set = dict_set(structure=structure, sort_structure=False)
+    incar_dict = {k.lower(): v for k, v in pmg_input_set.incar.items()}
+
+    potcar_symbols = pmg_input_set.potcar_symbols
+    potcar_setups = {symbol.split("_")[0]: symbol for symbol in potcar_symbols}
+    for k, v in potcar_setups.items():
+        if k in v:
+            potcar_setups[k] = v.split(k)[-1]
+
+    pp = pmg_input_set.potcar_functional.split("_")[0]
+
+    full_input_params = incar_dict | {"setups": potcar_setups, "pp": pp}
+
+    pmg_kpts = pmg_input_set.kpoints
+    if pmg_kpts is not None:
+        kpoints_dict = pmg_input_set.kpoints.as_dict()
+        full_input_params |= {
+            "kpts": kpoints_dict["kpoints"][0],
+            "gamma": kpoints_dict["generation_style"] == "Gamma",
+        }
+
+    return full_input_params, AseAtomsAdaptor().get_atoms(
+        pmg_input_set.poscar.structure
+    )
diff --git a/src/quacc/calculators/vasp/presets/MPScanSet.yaml b/src/quacc/calculators/vasp/presets/MPScanSet.yaml
deleted file mode 100644
index b2aaee4b..00000000
--- a/src/quacc/calculators/vasp/presets/MPScanSet.yaml
+++ /dev/null
@@ -1,30 +0,0 @@
-# Default VASP settings for calculations in the Materials Project
-# using the r2SCAN functional.
-# Caution: The pseudopotentials should be used with the .54 PAW PBE set
-# Source: https://github.com/materialsproject/atomate2/blob/main/src/atomate2/vasp/sets/BaseMPR2SCANRelaxSet.yaml
-# Last modified: 08-07-2023
-inputs:
-  algo: all
-  ediff: 1.0e-5
-  ediffg: -0.02
-  enaug: 1360
-  encut: 680
-  ibrion: 2
-  isif: 3
-  ismear: 0
-  kspacing: 0.22
-  laechg: True
-  lasph: True
-  lcharg: True
-  lmixtau: True
-  lorbit: 11
-  lreal: auto
-  lvtot: True
-  lwave: False
-  nelm: 200
-  nsw: 99
-  prec: accurate
-  sigma: 0.05
-  xc: r2scan
-parent_magmoms: magmoms_MP.yaml
-parent_setups: setups_pbe54_MP.yaml
diff --git a/src/quacc/calculators/vasp/presets/setups_pbe54_MP.yaml b/src/quacc/calculators/vasp/presets/setups_pbe54_MP.yaml
deleted file mode 100644
index 7d1b0ecc..00000000
--- a/src/quacc/calculators/vasp/presets/setups_pbe54_MP.yaml
+++ /dev/null
@@ -1,103 +0,0 @@
-# Materials Project pseudopotentials (r2SCAN workflow)
-# Note: These are specifically the .54 PAW PBE potentials
-# Source: https://github.com/materialsproject/atomate2/blob/main/src/atomate2/vasp/sets/BaseMPR2SCANRelaxSet.yaml
-# Last modified: 08-01-2023
-inputs:
-  setups:
-    Ac: Ac
-    Ag: Ag
-    Al: Al
-    Am: Am
-    Ar: Ar
-    As: As
-    At: At
-    Au: Au
-    B: B
-    Ba: Ba_sv
-    Be: Be_sv
-    Bi: Bi
-    Br: Br
-    C: C
-    Ca: Ca_sv
-    Cd: Cd
-    Ce: Ce
-    Cf: Cf
-    Cl: Cl
-    Cm: Cm
-    Co: Co
-    Cr: Cr_pv
-    Cs: Cs_sv
-    Cu: Cu_pv
-    Dy: Dy_3
-    Er: Er_3
-    Eu: Eu
-    F: F
-    Fe: Fe_pv
-    Fr: Fr_sv
-    Ga: Ga_d
-    Gd: Gd
-    Ge: Ge_d
-    H: H
-    He: He
-    Hf: Hf_pv
-    Hg: Hg
-    Ho: Ho_3
-    I: I
-    In: In_d
-    Ir: Ir
-    K: K_sv
-    Kr: Kr
-    La: La
-    Li: Li_sv
-    Lu: Lu_3
-    Mg: Mg_pv
-    Mn: Mn_pv
-    Mo: Mo_pv
-    N: N
-    Na: Na_pv
-    Nb: Nb_pv
-    Nd: Nd_3
-    Ne: Ne
-    Ni: Ni_pv
-    Np: Np
-    O: O
-    Os: Os_pv
-    P: P
-    Pa: Pa
-    Pb: Pb_d
-    Pd: Pd
-    Pm: Pm_3
-    Po: Po_d
-    Pr: Pr_3
-    Pt: Pt
-    Pu: Pu
-    Ra: Ra_sv
-    Rb: Rb_sv
-    Re: Re_pv
-    Rh: Rh_pv
-    Rn: Rn
-    Ru: Ru_pv
-    S: S
-    Sb: Sb
-    Sc: Sc_sv
-    Se: Se
-    Si: Si
-    Sm: Sm_3
-    Sn: Sn_d
-    Sr: Sr_sv
-    Ta: Ta_pv
-    Tb: Tb_3
-    Tc: Tc_pv
-    Te: Te
-    Th: Th
-    Ti: Ti_pv
-    Tl: Tl_d
-    Tm: Tm_3
-    U: U
-    V: V_pv
-    W: W_sv
-    Xe: Xe
-    Y: Y_sv
-    Yb: Yb_3
-    Zn: Zn
-    Zr: Zr_sv
diff --git a/src/quacc/calculators/vasp/vasp.py b/src/quacc/calculators/vasp/vasp.py
index d2c093b6..6861fe67 100644
--- a/src/quacc/calculators/vasp/vasp.py
+++ b/src/quacc/calculators/vasp/vasp.py
@@ -17,16 +17,19 @@ from quacc.calculators.vasp import vasp_custodian
 from quacc.calculators.vasp.io import load_vasp_yaml_calc
 from quacc.calculators.vasp.params import (
     get_param_swaps,
+    get_pmg_input_set_params,
     remove_unused_flags,
     set_auto_dipole,
     set_pmg_kpts,
 )
 from quacc.schemas.prep import set_magmoms
+from quacc.utils.dicts import sort_dict
 
 if TYPE_CHECKING:
     from typing import Literal
 
     from ase.atoms import Atoms
+    from pymatgen.io.vasp.sets import DictSet
 
 
 class Vasp(Vasp_):
@@ -50,6 +53,7 @@ class Vasp(Vasp_):
             | None
         ) = None,
         auto_dipole: bool | None = None,
+        pmg_input_set: DictSet | None = None,
         **kwargs,
     ) -> None:
         """
@@ -94,6 +98,9 @@ class Vasp(Vasp_):
         auto_dipole
             If True, will automatically set dipole moment correction parameters
             based on the center of mass (in the c dimension by default).
+        pmg_input_set
+            A Pymatgen input set to use for the VASP calculation, taken from a
+            `pymatgen.io.vasp.sets.DictSet` object.
         **kwargs
             Additional arguments to be passed to the VASP calculator, e.g.
             `xc='PBE'`, `encut=520`. Takes all valid ASE calculator arguments.
@@ -132,6 +139,7 @@ class Vasp(Vasp_):
         self.elemental_magmoms = elemental_magmoms
         self.pmg_kpts = pmg_kpts
         self.auto_dipole = auto_dipole
+        self.pmg_input_set = pmg_input_set
         self.kwargs = kwargs
 
         # Initialize for later
@@ -208,6 +216,14 @@ class Vasp(Vasp_):
             msg = "Atoms object has a constraint that is not compatible with Custodian."
             raise ValueError(msg)
 
+        # Get Pymatgen VASP input set parameters
+        if self.pmg_input_set:
+            pmg_calc_params, self.input_atoms = get_pmg_input_set_params(
+                self.pmg_input_set, self.input_atoms
+            )
+        else:
+            pmg_calc_params = {}
+
         # Get user-defined preset parameters for the calculator
         if self.preset:
             calc_preset = load_vasp_yaml_calc(SETTINGS.VASP_PRESET_DIR / self.preset)[
@@ -218,7 +234,7 @@ class Vasp(Vasp_):
 
         # Collect all the calculator parameters and prioritize the kwargs in the
         # case of duplicates.
-        self.user_calc_params = calc_preset | self.kwargs
+        self.user_calc_params = pmg_calc_params | calc_preset | self.kwargs
 
         # Allow the user to use setups='mysetups.yaml' to load in a custom
         # setups from a YAML file
@@ -271,4 +287,4 @@ class Vasp(Vasp_):
         )
 
         # Remove unused INCAR flags
-        self.user_calc_params = remove_unused_flags(self.user_calc_params)
+        self.user_calc_params = sort_dict(remove_unused_flags(self.user_calc_params))
diff --git a/src/quacc/recipes/vasp/mp.py b/src/quacc/recipes/vasp/mp.py
index 6d2b87eb..0178ef41 100644
--- a/src/quacc/recipes/vasp/mp.py
+++ b/src/quacc/recipes/vasp/mp.py
@@ -2,7 +2,6 @@
 Materials Project-compatible recipes.
 
 This set of recipes is meant to be compatible with the Materials Project
-Reference: https://doi.org/10.1103/PhysRevMaterials.6.013801
 
 !!! Info
 
@@ -16,10 +15,11 @@ Reference: https://doi.org/10.1103/PhysRevMaterials.6.013801
 
 from __future__ import annotations
 
+from functools import partial
 from pathlib import Path
 from typing import TYPE_CHECKING
 
-import numpy as np
+from pymatgen.io.vasp.sets import MPScanRelaxSet
 
 from quacc import flow, job
 from quacc.recipes.vasp._base import base_fn
@@ -30,14 +30,13 @@ if TYPE_CHECKING:
 
     from ase.atoms import Atoms
 
-    from quacc.schemas._aliases.vasp import MPRelaxFlowSchema, VaspSchema
+    from quacc.schemas._aliases.vasp import MPMetaGGARelaxFlowSchema, VaspSchema
 
 
 @job
-def mp_prerelax_job(
+def mp_metagga_prerelax_job(
     atoms: Atoms,
-    preset: str | None = "MPScanSet",
-    bandgap: float | None = None,
+    bandgap: float = 0.0,
     copy_files: str | Path | list[str | Path] | None = None,
     **calc_kwargs,
 ) -> VaspSchema:
@@ -45,12 +44,12 @@ def mp_prerelax_job(
     Function to pre-relax a structure with Materials Project settings. By default, this
     uses a PBEsol pre-relax step.
 
+    Reference: https://doi.org/10.1103/PhysRevMaterials.6.013801
+
     Parameters
     ----------
     atoms
         Atoms object
-    preset
-        Preset to use from `quacc.calculators.vasp.presets`.
     bandgap
         Estimate for the bandgap in eV.
     copy_files
@@ -68,27 +67,27 @@ def mp_prerelax_job(
     """
 
     calc_defaults = {
+        "pmg_input_set": partial(MPScanRelaxSet, bandgap=bandgap, auto_ismear=False),
         "ediffg": -0.05,
-        "xc": "pbesol",
+        "gga": "PS",
+        "laechg": False,  # Deviation from MP (but logical)
+        "lvtot": False,  # Deviation from MP (but logical)
         "lwave": True,
-        "lcharg": True,
-    } | _get_bandgap_swaps(bandgap)
-
+        "metagga": None,
+    }
     return base_fn(
         atoms,
-        preset=preset,
         calc_defaults=calc_defaults,
         calc_swaps=calc_kwargs,
-        additional_fields={"name": "MP Pre-Relax"},
+        additional_fields={"name": "MP Meta-GGA Pre-Relax"},
         copy_files=copy_files,
     )
 
 
 @job
-def mp_relax_job(
+def mp_metagga_relax_job(
     atoms: Atoms,
-    preset: str | None = "MPScanSet",
-    bandgap: float | None = None,
+    bandgap: float = 0.0,
     copy_files: str | Path | list[str | Path] | None = None,
     **calc_kwargs,
 ) -> VaspSchema:
@@ -96,12 +95,12 @@ def mp_relax_job(
     Function to relax a structure with Materials Project settings. By default, this uses
     an r2SCAN relax step.
 
+    Reference: https://doi.org/10.1103/PhysRevMaterials.6.013801
+
     Parameters
     ----------
     atoms
         Atoms object
-    preset
-        Preset to use from `quacc.calculators.vasp.presets`.
     bandgap
         Estimate for the bandgap in eV.
     copy_files
@@ -118,33 +117,91 @@ def mp_relax_job(
         See the type-hint for the data structure.
     """
 
-    calc_defaults = {"lcharg": True, "lwave": True} | _get_bandgap_swaps(bandgap)
+    calc_defaults = {
+        "pmg_input_set": partial(MPScanRelaxSet, bandgap=bandgap, auto_ismear=False),
+        "laechg": False,  # Deviation from MP (but logical)
+        "lvtot": False,  # Deviation from MP (but logical)
+        "lwave": True,
+    }
     return base_fn(
         atoms,
-        preset=preset,
         calc_defaults=calc_defaults,
         calc_swaps=calc_kwargs,
-        additional_fields={"name": "MP Relax"},
+        additional_fields={"name": "MP Meta-GGA Relax"},
+        copy_files=copy_files,
+    )
+
+
+@job
+def mp_metagga_static_job(
+    atoms: Atoms,
+    bandgap: float = 0.0,
+    copy_files: str | Path | list[str | Path] | None = None,
+    **calc_kwargs,
+) -> VaspSchema:
+    """
+    Function to run a static calculation on a structure with Materials Project settings. By default, this uses
+    an r2SCAN static step.
+
+    Parameters
+    ----------
+    atoms
+        Atoms object
+    bandgap
+        Estimate for the bandgap in eV.
+    copy_files
+        File(s) to copy to the runtime directory. If a directory is provided, it will be recursively unpacked.
+    **calc_kwargs
+        Dictionary of custom kwargs for the Vasp calculator. Set a value to
+        `None` to remove a pre-existing key entirely. For a list of available
+        keys, refer to `ase.calculators.vasp.vasp.Vasp`.
+
+    Returns
+    -------
+    VaspSchema
+        Dictionary of results from [quacc.schemas.vasp.vasp_summarize_run][].
+        See the type-hint for the data structure.
+    """
+
+    calc_defaults = {
+        "pmg_input_set": partial(MPScanRelaxSet, bandgap=bandgap, auto_ismear=False),
+        "algo": "fast",
+        "ismear": -5,
+        "lreal": False,
+        "lwave": True,  # Deviation from MP (but logical)
+        "nsw": 0,
+    }
+    return base_fn(
+        atoms,
+        calc_defaults=calc_defaults,
+        calc_swaps=calc_kwargs,
+        additional_fields={"name": "MP Meta-GGA Static"},
         copy_files=copy_files,
     )
 
 
 @flow
-def mp_relax_flow(
+def mp_metagga_relax_flow(
     atoms: Atoms,
     job_params: dict[str, dict[str, Any]] | None = None,
     job_decorators: dict[str, Callable | None] | None = None,
-) -> MPRelaxFlowSchema:
+) -> MPMetaGGARelaxFlowSchema:
     """
     Workflow consisting of:
 
     1. MP-compatible pre-relax
-        - name: "mp_prerelax_job"
-        - job: [quacc.recipes.vasp.mp.mp_prerelax_job][]
+        - name: "mp_metagga_prerelax_job"
+        - job: [quacc.recipes.vasp.mp.mp_metagga_prerelax_job][]
 
     2. MP-compatible relax
-        - name: "mp_relax_job"
-        - job: [quacc.recipes.vasp.mp.mp_relax_job][]
+        - name: "mp_metagga_relax_job"
+        - job: [quacc.recipes.vasp.mp.mp_metagga_relax_job][]
+
+    3. MP-compatible static
+        - name: "mp_metagga_static_job"
+        - job: [quacc.recipes.vasp.mp.mp_metagga_static_job][]
+
+    Reference: https://doi.org/10.1103/PhysRevMaterials.6.013801
 
     Parameters
     ----------
@@ -159,21 +216,25 @@ def mp_relax_flow(
 
     Returns
     -------
-    MPRelaxFlowSchema
+    MPMetaGGARelaxFlowSchema
         Dictionary of results. See the type-hint for the data structure.
     """
-    mp_prerelax_job_, mp_relax_job_ = customize_funcs(
-        ["mp_prerelax_job", "mp_relax_job"],
-        [mp_prerelax_job, mp_relax_job],
+    (
+        mp_metagga_prerelax_job_,
+        mp_metagga_relax_job_,
+        mp_metagga_static_job_,
+    ) = customize_funcs(
+        ["mp_metagga_prerelax_job", "mp_metagga_relax_job", "mp_metagga_static_job"],
+        [mp_metagga_prerelax_job, mp_metagga_relax_job, mp_metagga_static_job],
         parameters=job_params,
         decorators=job_decorators,
     )
 
     # Run the prerelax
-    prerelax_results = mp_prerelax_job_(atoms)
+    prerelax_results = mp_metagga_prerelax_job_(atoms)
 
     # Run the relax
-    relax_results = mp_relax_job_(
+    relax_results = mp_metagga_relax_job_(
         prerelax_results["atoms"],
         bandgap=prerelax_results["output"]["bandgap"],
         copy_files=[
@@ -183,28 +244,18 @@ def mp_relax_flow(
     )
     relax_results["prerelax"] = prerelax_results
 
-    return relax_results
-
-
-def _get_bandgap_swaps(bandgap: float | None = None) -> dict[str, float]:
-    """
-    Get bandgap-related swaps.
-
-    Parameters
-    ---------
-    bandgap
-        The bandgap, in units of eV.
-
-    Returns
-    -------
-    dict
-        Dictionary of swaps.
-    """
+    # Run the static
+    static_results = mp_metagga_static_job_(
+        relax_results["atoms"],
+        bandgap=relax_results["output"]["bandgap"],
+        copy_files=[
+            Path(relax_results["dir_name"]) / "CHGCAR",
+            Path(relax_results["dir_name"]) / "WAVECAR",
+        ],
+    )
 
-    if bandgap is None:
-        return {"kspacing": 0.22, "ismear": 0, "sigma": 0.05}
-    if bandgap <= 1e-4:
-        return {"kspacing": 0.22, "ismear": 2, "sigma": 0.2}
-    rmin = max(1.5, 25.22 - 2.87 * bandgap)
-    kspacing = 2 * np.pi * 1.0265 / (rmin - 1.0183)
-    return {"kspacing": min(kspacing, 0.44), "ismear": -5, "sigma": 0.05}
+    return {
+        "prerelax": prerelax_results,
+        "relax": relax_results,
+        "static": static_results,
+    }
diff --git a/src/quacc/schemas/_aliases/vasp.py b/src/quacc/schemas/_aliases/vasp.py
index 5db65d7d..00634f2d 100644
--- a/src/quacc/schemas/_aliases/vasp.py
+++ b/src/quacc/schemas/_aliases/vasp.py
@@ -58,10 +58,12 @@ class DoubleRelaxSchema(VaspSchema):
     relax1: VaspSchema
 
 
-class MPRelaxFlowSchema(VaspSchema):
-    """Type hint associated with the MP relaxation flows."""
+class MPMetaGGARelaxFlowSchema(VaspSchema):
+    """Type hint associated with the MP meta-GGA relaxation flows."""
 
     prerelax: VaspSchema
+    relax: VaspSchema
+    static: VaspSchema
 
 
 class QMOFRelaxSchema(VaspSchema):


Test patch:
diff --git a/tests/core/calculators/vasp/test_vasp.py b/tests/core/calculators/vasp/test_vasp.py
index 9f35f07a..1069e6e0 100644
--- a/tests/core/calculators/vasp/test_vasp.py
+++ b/tests/core/calculators/vasp/test_vasp.py
@@ -9,6 +9,7 @@ from ase.calculators.singlepoint import SinglePointDFTCalculator
 from ase.calculators.vasp import Vasp as Vasp_
 from ase.constraints import FixAtoms, FixBondLength
 from ase.io import read
+from pymatgen.io.vasp.sets import MPRelaxSet, MPScanRelaxSet
 
 from quacc import SETTINGS
 from quacc.calculators.vasp import Vasp, presets
@@ -69,9 +70,9 @@ def test_presets():
     assert calc.exp_params["ediff"] == 1e-5
     assert calc.float_params["encut"] == 450
 
-    calc = Vasp(atoms, xc="scan", preset="MPScanSet")
+    calc = Vasp(atoms, xc="scan", pmg_input_set=MPScanRelaxSet)
     assert calc.xc.lower() == "scan"
-    assert calc.string_params["algo"] == "all"
+    assert calc.string_params["algo"].lower() == "all"
     assert calc.exp_params["ediff"] == 1e-5
 
 
@@ -164,6 +165,7 @@ def test_magmoms(atoms_mag, atoms_nomag, atoms_nospin):
     atoms[-1].symbol = "Fe"
     calc = Vasp(atoms, preset="BulkSet")
     atoms.calc = calc
+    assert atoms.get_chemical_symbols() == ["Cu", "Cu", "Cu", "Fe"]
     assert atoms.get_initial_magnetic_moments().tolist() == [2.0] * (len(atoms) - 1) + [
         5.0
     ]
@@ -186,16 +188,15 @@ def test_magmoms(atoms_mag, atoms_nomag, atoms_nospin):
 
     atoms = bulk("Cu") * (2, 2, 1)
     atoms[-1].symbol = "Fe"
-    calc = Vasp(atoms, preset="MPScanSet")
+    calc = Vasp(atoms, pmg_input_set=MPScanRelaxSet)
     atoms.calc = calc
-    assert atoms.get_initial_magnetic_moments().tolist() == [1.0] * (len(atoms) - 1) + [
-        5.0
-    ]
+    assert atoms.get_chemical_symbols() == ["Cu", "Cu", "Cu", "Fe"]
+    assert calc.parameters["magmom"] == [0.6, 0.6, 0.6, 5.0]
 
     atoms = bulk("Cu") * (2, 2, 1)
     atoms[-1].symbol = "Fe"
     atoms.set_initial_magnetic_moments([3.14] * (len(atoms) - 1) + [1.0])
-    calc = Vasp(atoms, preset="BulkSet")
+    calc = Vasp(atoms, pmg_input_set=MPScanRelaxSet)
     atoms.calc = calc
     assert atoms.get_initial_magnetic_moments().tolist() == [3.14] * (
         len(atoms) - 1
@@ -707,7 +708,7 @@ def test_setups():
     assert calc.parameters["setups"]["Cu"] == ""
 
     atoms = bulk("Cu")
-    calc = Vasp(atoms, preset="MPScanSet")
+    calc = Vasp(atoms, pmg_input_set=MPScanRelaxSet)
     assert calc.parameters["setups"]["Cu"] == "_pv"
 
     atoms = bulk("Cu")
@@ -719,7 +720,12 @@ def test_setups():
     assert calc.parameters["setups"]["Cu"] == ""
 
     atoms = bulk("Cu")
-    calc = Vasp(atoms, setups="minimal", preset="MPScanSet")
+    calc = Vasp(atoms, setups="minimal", preset="BulkSet")
+    assert isinstance(calc.parameters["setups"], str)
+    assert calc.parameters["setups"] == "minimal"
+
+    atoms = bulk("Cu")
+    calc = Vasp(atoms, setups="minimal", pmg_input_set=MPScanRelaxSet)
     assert isinstance(calc.parameters["setups"], str)
     assert calc.parameters["setups"] == "minimal"
 
@@ -823,3 +829,66 @@ def test_preset_override():
 
     calc = Vasp(atoms, preset="BulkSet", efermi=None)
     assert calc.parameters.get("efermi") is None
+
+
+def test_pmg_input_set():
+    atoms = bulk("Cu")
+    calc = Vasp(atoms, pmg_input_set=MPRelaxSet, incar_copilot="off")
+    assert calc.parameters == {
+        "algo": "Fast",
+        "ediff": 5e-05,
+        "encut": 520,
+        "ibrion": 2,
+        "isif": 3,
+        "ismear": -5,
+        "ispin": 2,
+        "lasph": True,
+        "lorbit": 11,
+        "lreal": "Auto",
+        "lwave": False,
+        "nelm": 100,
+        "nsw": 99,
+        "pp": "PBE",
+        "prec": "Accurate",
+        "sigma": 0.05,
+        "magmom": [0.6],
+        "lmaxmix": 4,
+        "kpts": [11, 11, 11],
+        "gamma": True,
+        "setups": {"Cu": "_pv"},
+    }
+
+
+def test_pmg_input_set2():
+    atoms = bulk("Fe") * (2, 1, 1)
+    atoms[0].symbol = "O"
+    calc = Vasp(atoms, pmg_input_set=MPRelaxSet, incar_copilot="off")
+    assert calc.parameters == {
+        "algo": "Fast",
+        "ediff": 0.0001,
+        "encut": 520,
+        "ibrion": 2,
+        "isif": 3,
+        "ismear": -5,
+        "ispin": 2,
+        "lasph": True,
+        "ldau": True,
+        "ldauj": [0, 0],
+        "ldaul": [0, 2.0],
+        "ldautype": 2,
+        "ldauu": [0, 5.3],
+        "ldauprint": 1,
+        "lorbit": 11,
+        "lreal": "Auto",
+        "lwave": False,
+        "nelm": 100,
+        "nsw": 99,
+        "pp": "PBE",
+        "prec": "Accurate",
+        "sigma": 0.05,
+        "magmom": [2.3, 2.3],
+        "lmaxmix": 4,
+        "kpts": [5, 11, 11],
+        "gamma": True,
+        "setups": {"Fe": "_pv", "O": ""},
+    }
diff --git a/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py b/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py
index 035a83e7..65e433fa 100644
--- a/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py
+++ b/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py
@@ -3,7 +3,12 @@ from ase.build import bulk, molecule
 
 from quacc import SETTINGS
 from quacc.recipes.vasp.core import double_relax_job, relax_job, static_job
-from quacc.recipes.vasp.mp import mp_prerelax_job, mp_relax_flow, mp_relax_job
+from quacc.recipes.vasp.mp import (
+    mp_metagga_prerelax_job,
+    mp_metagga_relax_flow,
+    mp_metagga_relax_job,
+    mp_metagga_static_job,
+)
 from quacc.recipes.vasp.qmof import qmof_relax_job
 from quacc.recipes.vasp.slabs import bulk_to_slabs_flow
 from quacc.recipes.vasp.slabs import relax_job as slab_relax_job
@@ -351,105 +356,208 @@ def test_qmof(tmp_path, monkeypatch):
     assert output["double_relax"][1]["parameters"]["isif"] == 2
 
 
-def test_mp_prerelax_job(tmp_path, monkeypatch):
+def test_mp_metagga_prerelax_job(tmp_path, monkeypatch):
     monkeypatch.chdir(tmp_path)
 
     atoms = bulk("Al")
-    output = mp_prerelax_job(atoms)
+    output = mp_metagga_prerelax_job(atoms)
     assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "pbesol"
+    assert output["parameters"] == {
+        "algo": "All",
+        "ediff": 1e-5,
+        "ediffg": -0.05,
+        "efermi": "midgap",  # added by copilot
+        "enaug": 1360,
+        "encut": 680,
+        "gga": "PS",
+        "ibrion": 2,
+        "isif": 3,
+        "ismear": 0,
+        "ispin": 2,
+        "kspacing": 0.22,
+        "laechg": False,  # disabled by us
+        "lasph": True,
+        "lcharg": True,
+        "lelf": False,
+        "lmixtau": True,
+        "lorbit": 11,
+        "lreal": "Auto",
+        "lvtot": False,  # disabled by us
+        "lwave": True,
+        "magmom": [0.6],
+        "nelm": 200,
+        "nsw": 99,
+        "prec": "Accurate",
+        "setups": {"Al": ""},
+        "sigma": 0.05,
+        "pp": "PBE",
+    }
+
+    output = mp_metagga_prerelax_job(atoms, bandgap=0)
+    assert output["nsites"] == len(atoms)
+    assert output["parameters"]["gga"] == "PS"
     assert output["parameters"]["ediffg"] == -0.05
     assert output["parameters"]["encut"] == 680
     assert output["parameters"]["kspacing"] == 0.22
     assert output["parameters"]["ismear"] == 0
     assert output["parameters"]["sigma"] == 0.05
+    assert output["parameters"]["pp"] == "PBE"
+    assert "metagga" not in output["parameters"]
 
-    output = mp_prerelax_job(atoms, bandgap=0)
-    assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "pbesol"
-    assert output["parameters"]["ediffg"] == -0.05
-    assert output["parameters"]["encut"] == 680
-    assert output["parameters"]["kspacing"] == 0.22
-    assert output["parameters"]["ismear"] == 2
-    assert output["parameters"]["sigma"] == 0.2
-
-    output = mp_prerelax_job(atoms, bandgap=100)
+    output = mp_metagga_prerelax_job(atoms, bandgap=100)
     assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "pbesol"
+    assert output["parameters"]["gga"] == "PS"
     assert output["parameters"]["ediffg"] == -0.05
     assert output["parameters"]["encut"] == 680
     assert output["parameters"]["kspacing"] == 0.44
-    assert output["parameters"]["ismear"] == -5
+    assert output["parameters"]["ismear"] == 0
     assert output["parameters"]["sigma"] == 0.05
+    assert output["parameters"]["pp"] == "PBE"
+    assert "metagga" not in output["parameters"]
 
 
-def test_mp_relax_job(tmp_path, monkeypatch):
+def test_mp_metagga_relax_job(tmp_path, monkeypatch):
     monkeypatch.chdir(tmp_path)
 
     atoms = bulk("Al")
 
-    output = mp_relax_job(atoms)
+    output = mp_metagga_relax_job(atoms)
+    assert output["nsites"] == len(atoms)
+    assert output["parameters"] == {
+        "algo": "All",
+        "ediff": 1e-5,
+        "ediffg": -0.02,
+        "efermi": "midgap",  # added by copilot
+        "enaug": 1360,
+        "encut": 680,
+        "ibrion": 2,
+        "isif": 3,
+        "ismear": 0,
+        "ispin": 2,
+        "kspacing": 0.22,
+        "laechg": False,  # disabled by us
+        "lasph": True,
+        "lcharg": True,
+        "lelf": False,
+        "lmixtau": True,
+        "lorbit": 11,
+        "lreal": "Auto",
+        "lvtot": False,  # disabled by us
+        "lwave": True,
+        "magmom": [0.6],
+        "metagga": "R2scan",
+        "nelm": 200,
+        "nsw": 99,
+        "prec": "Accurate",
+        "sigma": 0.05,
+        "pp": "PBE",
+        "setups": {"Al": ""},
+    }
+
+    output = mp_metagga_relax_job(atoms, bandgap=0)
     assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "r2scan"
+    assert output["parameters"]["metagga"].lower() == "r2scan"
     assert output["parameters"]["ediffg"] == -0.02
     assert output["parameters"]["encut"] == 680
     assert output["parameters"]["kspacing"] == 0.22
     assert output["parameters"]["ismear"] == 0
     assert output["parameters"]["sigma"] == 0.05
+    assert output["parameters"]["pp"] == "PBE"
 
-    output = mp_relax_job(atoms, bandgap=0)
+    output = mp_metagga_relax_job(atoms, bandgap=100)
     assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "r2scan"
-    assert output["parameters"]["ediffg"] == -0.02
-    assert output["parameters"]["encut"] == 680
-    assert output["parameters"]["kspacing"] == 0.22
-    assert output["parameters"]["ismear"] == 2
-    assert output["parameters"]["sigma"] == 0.2
-
-    output = mp_relax_job(atoms, bandgap=100)
-    assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "r2scan"
+    assert output["parameters"]["metagga"].lower() == "r2scan"
     assert output["parameters"]["ediffg"] == -0.02
     assert output["parameters"]["encut"] == 680
     assert output["parameters"]["kspacing"] == 0.44
-    assert output["parameters"]["ismear"] == -5
+    assert output["parameters"]["ismear"] == 0
     assert output["parameters"]["sigma"] == 0.05
+    assert output["parameters"]["pp"] == "PBE"
 
 
-def test_mp_relax_flow(tmp_path, monkeypatch):
+def test_mp_metagga_static_job(tmp_path, monkeypatch):
     monkeypatch.chdir(tmp_path)
 
     atoms = bulk("Al")
 
-    output = mp_relax_flow(atoms)
+    output = mp_metagga_static_job(atoms)
     assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "r2scan"
-    assert output["parameters"]["ediffg"] == -0.02
-    assert output["parameters"]["encut"] == 680
-    assert output["parameters"]["ismear"] == 2
-    assert output["parameters"]["sigma"] == 0.2
-    assert output["parameters"]["kspacing"] == 0.22
-    assert output["prerelax"]["parameters"]["xc"] == "pbesol"
+    assert output["parameters"] == {
+        "algo": "fast",
+        "ediff": 1e-05,
+        "efermi": "midgap",  # added by copilot
+        "enaug": 1360,
+        "encut": 680,
+        "ismear": -5,
+        "ispin": 2,
+        "kspacing": 0.22,
+        "laechg": True,
+        "lasph": True,
+        "lcharg": True,
+        "lelf": False,
+        "lmixtau": True,
+        "lorbit": 11,
+        "lreal": False,
+        "lvtot": True,
+        "lwave": True,  # enabled by us
+        "magmom": [0.6],
+        "metagga": "R2scan",
+        "nelm": 200,
+        "nsw": 0,
+        "prec": "Accurate",
+        "sigma": 0.05,
+        "pp": "PBE",
+        "setups": {"Al": ""},
+    }
+
+
+def test_mp_metagga_relax_flow(tmp_path, monkeypatch):
+    monkeypatch.chdir(tmp_path)
+
+    atoms = bulk("Al")
+
+    output = mp_metagga_relax_flow(atoms)
+    assert output["static"]["nsites"] == len(atoms)
+    assert output["relax"]["parameters"]["metagga"].lower() == "r2scan"
+    assert output["relax"]["parameters"]["ediffg"] == -0.02
+    assert output["relax"]["parameters"]["encut"] == 680
+    assert output["relax"]["parameters"]["ismear"] == 0
+    assert output["relax"]["parameters"]["sigma"] == 0.05
+    assert output["relax"]["parameters"]["kspacing"] == 0.22
+    assert output["relax"]["parameters"]["pp"] == "PBE"
+    assert output["prerelax"]["parameters"]["gga"] == "PS"
     assert output["prerelax"]["parameters"]["ismear"] == 0
+    assert output["prerelax"]["parameters"]["pp"] == "PBE"
 
     atoms = bulk("C")
-    output = mp_relax_flow(atoms)
-    assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "r2scan"
-    assert output["parameters"]["ediffg"] == -0.02
-    assert output["parameters"]["encut"] == 680
-    assert output["parameters"]["ismear"] == -5
-    assert output["parameters"]["kspacing"] == pytest.approx(0.28329488761304206)
+    output = mp_metagga_relax_flow(atoms)
+    assert output["static"]["nsites"] == len(atoms)
+    assert output["relax"]["parameters"]["metagga"].lower() == "r2scan"
+    assert output["relax"]["parameters"]["ediffg"] == -0.02
+    assert output["relax"]["parameters"]["encut"] == 680
+    assert output["relax"]["parameters"]["ismear"] == 0
+    assert output["relax"]["parameters"]["kspacing"] == pytest.approx(
+        0.28329488761304206
+    )
+    assert output["relax"]["parameters"]["pp"] == "PBE"
     assert output["prerelax"]["parameters"]["ismear"] == 0
+    assert output["prerelax"]["parameters"]["pp"] == "PBE"
 
     atoms = molecule("O2")
     atoms.center(vacuum=10)
     atoms.pbc = True
-    output = mp_relax_flow(atoms)
-    assert output["nsites"] == len(atoms)
-    assert output["parameters"]["xc"] == "r2scan"
-    assert output["parameters"]["ediffg"] == -0.02
-    assert output["parameters"]["encut"] == 680
-    assert output["parameters"]["ismear"] == -5
-    assert output["parameters"]["kspacing"] == pytest.approx(0.28329488761304206)
+    output = mp_metagga_relax_flow(atoms)
+    assert output["static"]["nsites"] == len(atoms)
+    assert output["static"]["parameters"]["ismear"] == -5
+    assert output["static"]["parameters"]["nsw"] == 0
+    assert output["static"]["parameters"]["algo"] == "fast"
+    assert output["relax"]["parameters"]["metagga"].lower() == "r2scan"
+    assert output["relax"]["parameters"]["ediffg"] == -0.02
+    assert output["relax"]["parameters"]["encut"] == 680
+    assert output["relax"]["parameters"]["ismear"] == 0
+    assert output["relax"]["parameters"]["kspacing"] == pytest.approx(
+        0.28329488761304206
+    )
+    assert output["relax"]["parameters"]["pp"] == "PBE"
     assert output["prerelax"]["parameters"]["ismear"] == 0
+    assert output["prerelax"]["parameters"]["pp"] == "PBE"
diff --git a/tests/requirements.txt b/tests/requirements.txt
index ee74d6d0..cc5ff639 100644
--- a/tests/requirements.txt
+++ b/tests/requirements.txt
@@ -8,5 +8,5 @@ numpy==1.26.4
 psutil==5.9.8
 pydantic==2.6.1
 pydantic-settings==2.1.0
-pymatgen==2024.2.8
+git+https://github.com/materialsproject/pymatgen.git
 typer[all]==0.9.0