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.

This workflow was used to build the Materials Project’s Harmonic Phonon Database, described in this preprint.

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[pheasy]'. It includes phonopy and ALM. ALM is compiled from source. See the atomate2 VASP documentation if that build fails.

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, third order by default, set by anhar_max_order) and fcs_cutoff_radius (cutoff radii in Bohr for the third- and fourth-order force constants, the first entry is not used). 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 the displaced supercells (a single finite displacement for silicon), 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

To go beyond the harmonic approximation, set cal_anhar_fcs=True. Pheasy will then extract third-order force constants with LASSO regression, and fourth-order ones too with anhar_max_order=4 (the cutoff radii are controlled by fcs_cutoff_radius). The number of anharmonic displaced supercells is set from the number of free force constants, with at least 20. Above 600 the job stops unless num_disp_anhar is set.

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.