This first part is only needed as we have to mock VASP here as we cannot run it directly in a jupyter notebook:

import warnings

from mock_vasp import TEST_DIR, mock_vasp

ref_paths = {
    "tight relax 1": "Si_pheasy/tight_relax_1",
    "tight relax 2": "Si_pheasy/tight_relax_2",
    "phonon static 1/2": "Si_pheasy/phonon_static_1_2",
    "phonon static 2/2": "Si_pheasy/phonon_static_2_2",
    "static": "Si_pheasy/static",
    "dielectric": "Si_pheasy/dielectric",
}

Lattice-dynamics workflow with Pheasy and VASP

Background

The Pheasy code extracts (an)harmonic interatomic force constants from a set of displaced supercell calculations using machine-learning (LASSO) regression. Compared to the finite-displacement approach in the standard phonon workflow, far fewer (randomly displaced) supercells are needed, which substantially reduces the number of DFT calculations.

The workflow has the same basic structure as the phonopy-based phonon workflow and uses Phonopy to post-process the force constants into band structures, densities of states, and thermodynamic properties.

To run this tutorial, the pheasy extra must be installed: pip install 'atomate2[phonons,pheasy]'. Anharmonic force-constant extraction (cal_anhar_fcs=True) additionally requires the alamode extra — see the atomate2 VASP documentation for installation hints.

Let’s run the workflow

First we load the structure and everything needed to run the workflow.

from jobflow import JobStore, run_locally
from maggma.stores import MemoryStore
from pymatgen.core import Structure

from atomate2.vasp.flows.pheasy import PhononMaker
from atomate2.vasp.powerups import update_user_incar_settings

warnings.filterwarnings("ignore")

job_store = JobStore(MemoryStore(), additional_stores={"data": MemoryStore()})
si_structure = Structure.from_file(
    TEST_DIR / "vasp" / "Si_pheasy" / "tight_relax_1" / "inputs" / "POSCAR.gz"
)

Then we use the pheasy PhononMaker to generate a Flow. The maker mirrors the phonopy-based PhononMaker; the most important additional switches are cal_anhar_fcs (extract anharmonic force constants up to fourth order and renormalize the phonon energies) and fcs_cutoff_radius (cutoff radii for the second-, third- and fourth-order force constants). Here we only extract harmonic force constants for silicon. As always, a tight structural relaxation is performed first — make sure it is converged very accurately for production runs.

flow = PhononMaker(
    force_diagonal=True,
    min_length=12,
    cal_anhar_fcs=False,
    create_thermal_displacements=True,
).make(structure=si_structure)

flow = update_user_incar_settings(
    flow,
    {
        "ENCUT": 600,
        "ISMEAR": 0,
        "SIGMA": 0.05,
        "KSPACING": 0.15,
        "ISPIN": 1,
        "EDIFFG": -1e-04,
        "EDIFF": 1e-07,
    },
)

The flow relaxes the bulk structure, generates a set of randomly displaced supercells, computes their forces with VASP, and then fits the force constants with pheasy. We can visualize the flow first.

flow.draw_graph().show()

We now run the flow with run_locally. We mock the VASP runs here. Normally, you would simply use run_locally without the with mock_vasp context manager.

with mock_vasp(ref_paths=ref_paths) as mf:
    run_locally(
        flow,
        create_folders=True,
        ensure_success=True,
        raise_immediately=True,
        store=job_store,
    )

Let’s retrieve the phonon band structure and density of states from the job store and plot them.

from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter

job_store.connect()

result = job_store.query_one(
    {"name": "generate_frequencies_eigenvectors"},
    properties=[
        "output.phonon_dos",
        "output.phonon_bandstructure",
    ],
    load=True,
    sort={"completed_at": -1},  # to get the latest computation
)

from emmet.core.phonon import PhononBS, PhononDOS

ph_bs = PhononBS(
    **result["output"]["phonon_bandstructure"]
).to_pmg  # get pymatgen bandstructure object
ph_dos = PhononDOS(
    **result["output"]["phonon_dos"]
).to_pmg  # get pymatgen phonon dos object

# initialize dos plotter and visualize dos plot
dos_plot = PhononDosPlotter()
dos_plot.add_dos(label="a", dos=ph_dos)
dos_plot.get_plot()

# initialize Phonon bandstructure plotter and visualize band structure plot
bs_plot = PhononBSPlotter(bs=ph_bs)
bs_plot.get_plot()

Anharmonic force constants and phonon renormalization

To go beyond the harmonic approximation, set cal_anhar_fcs=True. Pheasy will then extract third- and fourth-order force constants with LASSO regression (the cutoff radii are controlled by fcs_cutoff_radius). With renorm_phonon=True, the phonon energies are additionally renormalized at the temperatures given by renorm_temp. The lattice thermal conductivity can additionally be computed with cal_ther_cond=True (requires phono3py). Both options need more displaced supercells than the harmonic run shown here — pheasy will tell you how many it generated.

The same workflow is also available for forcefields via from atomate2.forcefields.flows.pheasy import PhononMaker, which is a cheap way to try it out without DFT.