Skip to content

Commit

Permalink
Merge branch 'main' of github.com:materialsproject/emmet into main
Browse files Browse the repository at this point in the history
  • Loading branch information
@tschaume
tschaume committed Mar 24, 2023
2 parents c929714 + 0b291b5 commit 83b9b95
Showing 1 changed file with 157 additions and 37 deletions.
194 changes: 157 additions & 37 deletions emmet-builders/emmet/builders/abinit/phonon.py
View file
@@ -1,8 +1,10 @@
import tempfile
import os
from math import ceil
from emmet.builders.settings import EmmetBuildSettings
import numpy as np
from typing import Optional, Dict, List, Iterator, Tuple
from maggma.utils import grouper

from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine
from pymatgen.phonon.dos import CompletePhononDos
Expand All @@ -20,8 +22,19 @@
from maggma.builders import Builder
from maggma.core import Store

from emmet.core.phonon import PhononWarnings, ThermodynamicProperties, AbinitPhonon, VibrationalEnergy
from emmet.core.phonon import PhononDos, PhononBandStructure, PhononWebsiteBS, Ddb, ThermalDisplacement
from emmet.core.phonon import (
PhononWarnings,
ThermodynamicProperties,
AbinitPhonon,
VibrationalEnergy,
)
from emmet.core.phonon import (
PhononDos,
PhononBandStructure,
PhononWebsiteBS,
Ddb,
ThermalDisplacement,
)
from emmet.core.polar import DielectricDoc, BornEffectiveCharges, IRDielectric
from emmet.core.utils import jsanitize

Expand Down Expand Up @@ -105,6 +118,25 @@ def __init__(
**kwargs,
)

def prechunk(self, number_splits: int): # pragma: no cover
"""
Gets all materials that need phonons
Returns:
generator of materials to extract phonon properties
"""

# All relevant materials that have been updated since phonon props were last calculated
q = dict(self.query)

mats = self.phonon.newer_in(self.phonon_materials, exhaustive=True, criteria=q)

N = ceil(len(mats) / number_splits)

for mpid_chunk in grouper(mats, N):

yield {"query": {self.phonon.key: {"$in": list(mpid_chunk)}}}

def get_items(self) -> Iterator[Dict]:
"""
Gets all materials that need phonons
Expand Down Expand Up @@ -134,19 +166,27 @@ def get_items(self) -> Iterator[Dict]:
}

for m in mats:
item = self.phonon_materials.query_one(properties=projection, criteria={self.phonon_materials.key: m})
item = self.phonon_materials.query_one(
properties=projection, criteria={self.phonon_materials.key: m}
)

# Read the DDB file and pass as an object. Do not write here since in case of parallel
# execution each worker will write its own file.
ddb_data = self.ddb_source.query_one(criteria={"_id": item["abinit_output"]["ddb_id"]})
ddb_data = self.ddb_source.query_one(
criteria={"_id": item["abinit_output"]["ddb_id"]}
)
if not ddb_data:
self.logger.warning(f"DDB file not found for file id {item['abinit_output']['ddb_id']}")
self.logger.warning(
f"DDB file not found for file id {item['abinit_output']['ddb_id']}"
)
continue

try:
item["ddb_str"] = ddb_data["data"].decode("utf-8")
except Exception:
self.logger.warning(f"could not extract DDB for file id {item['abinit_output']['ddb_id']}")
self.logger.warning(
f"could not extract DDB for file id {item['abinit_output']['ddb_id']}"
)
continue

yield item
Expand All @@ -170,9 +210,13 @@ def process_item(self, item: Dict) -> Optional[Dict]:
abinit_input_vars = self.abinit_input_vars(item)
phonon_properties = self.get_phonon_properties(item)
sr_break = self.get_sum_rule_breakings(item)
ph_warnings = get_warnings(sr_break["asr"], sr_break["cnsr"], phonon_properties["ph_bs"])
ph_warnings = get_warnings(
sr_break["asr"], sr_break["cnsr"], phonon_properties["ph_bs"]
)
if PhononWarnings.NEG_FREQ not in ph_warnings:
thermodynamic, vibrational_energy = get_thermodynamic_properties(phonon_properties["ph_dos"])
thermodynamic, vibrational_energy = get_thermodynamic_properties(
phonon_properties["ph_dos"]
)
else:
thermodynamic, vibrational_energy = None, None

Expand Down Expand Up @@ -201,10 +245,14 @@ def process_item(self, item: Dict) -> Optional[Dict]:
abinit_input_vars=abinit_input_vars,
)

phbs = PhononBandStructure(material_id=item["mp_id"], band_structure=phonon_properties["ph_bs"].as_dict())
phbs = PhononBandStructure(
material_id=item["mp_id"],
band_structure=phonon_properties["ph_bs"].as_dict(),
)

phws = PhononWebsiteBS(
material_id=item["mp_id"], phononwebsite=phonon_properties["ph_bs"].as_phononwebsite()
material_id=item["mp_id"],
phononwebsite=phonon_properties["ph_bs"].as_phononwebsite(),
)

phdos = PhononDos(
Expand All @@ -213,7 +261,7 @@ def process_item(self, item: Dict) -> Optional[Dict]:
dos_method=phonon_properties["ph_dos_method"],
)

ddb = Ddb(material_id=item["mp_id"], ddb=item["ddb_str"],)
ddb = Ddb(material_id=item["mp_id"], ddb=item["ddb_str"])

th_disp = ThermalDisplacement(
material_id=item["mp_id"],
Expand Down Expand Up @@ -242,7 +290,11 @@ def process_item(self, item: Dict) -> Optional[Dict]:

return d
except Exception as error:
self.logger.warning("Error generating the phonon properties for {}: {}".format(item["mp_id"], error))
self.logger.warning(
"Error generating the phonon properties for {}: {}".format(
item["mp_id"], error
)
)
return None

def get_phonon_properties(self, item: Dict) -> Dict:
Expand All @@ -266,11 +318,19 @@ def get_phonon_properties(self, item: Dict) -> Dict:
has_epsinf = ddb.has_epsinf_terms()

anaddb_input, labels_list = self.get_properties_anaddb_input(
item, bs=True, dos="tetra", lo_to_splitting=has_bec, use_dieflag=has_epsinf
item,
bs=True,
dos="tetra",
lo_to_splitting=has_bec,
use_dieflag=has_epsinf,
)
task = self.run_anaddb(
ddb_path=ddb_path, anaddb_input=anaddb_input, workdir=workdir
)
task = self.run_anaddb(ddb_path=ddb_path, anaddb_input=anaddb_input, workdir=workdir,)

with task.open_phbst() as phbst_file, AnaddbNcFile(task.outpath_from_ext("anaddb.nc")) as ananc_file:
with task.open_phbst() as phbst_file, AnaddbNcFile(
task.outpath_from_ext("anaddb.nc")
) as ananc_file:
# phbst
phbands = phbst_file.phbands
if has_bec:
Expand All @@ -286,7 +346,9 @@ def get_phonon_properties(self, item: Dict) -> Dict:
e_electronic = ananc_file.epsinf.tolist()
else:
e_electronic = None
e_total = ananc_file.eps0.tolist() if ananc_file.eps0 is not None else None
e_total = (
ananc_file.eps0.tolist() if ananc_file.eps0 is not None else None
)
if e_electronic and e_total:
e_ionic = (ananc_file.eps0 - ananc_file.epsinf).tolist()
dielectric = DielectricDoc.from_ionic_and_electronic(
Expand All @@ -300,11 +362,20 @@ def get_phonon_properties(self, item: Dict) -> Dict:
dielectric = None

# both
if e_electronic and e_total and ananc_file.oscillator_strength is not None:
die_gen = DielectricTensorGenerator.from_objects(phbands, ananc_file)
if (
e_electronic
and e_total
and ananc_file.oscillator_strength is not None
):
die_gen = DielectricTensorGenerator.from_objects(
phbands, ananc_file
)

ir_tensor = IRDielectricTensor(
die_gen.oscillator_strength, die_gen.phfreqs, die_gen.epsinf, die_gen.structure
die_gen.oscillator_strength,
die_gen.phfreqs,
die_gen.epsinf,
die_gen.structure,
).as_dict()
ir_spectra = IRDielectric(ir_dielectric_tensor=ir_tensor)
else:
Expand All @@ -328,9 +399,17 @@ def get_phonon_properties(self, item: Dict) -> Dict:
)
with tempfile.TemporaryDirectory() as workdir_dos:
anaddb_input_dos, _ = self.get_properties_anaddb_input(
item, bs=False, dos="gauss", lo_to_splitting=has_bec, use_dieflag=has_epsinf
item,
bs=False,
dos="gauss",
lo_to_splitting=has_bec,
use_dieflag=has_epsinf,
)
task_dos = self.run_anaddb(
ddb_path=ddb_path,
anaddb_input=anaddb_input_dos,
workdir=workdir_dos,
)
task_dos = self.run_anaddb(ddb_path=ddb_path, anaddb_input=anaddb_input_dos, workdir=workdir_dos)
with task_dos.open_phdos() as phdos_file:
complete_dos = phdos_file.to_pymatgen()
msqd_dos = phdos_file.msqd_dos
Expand All @@ -354,7 +433,9 @@ def get_sum_rule_breakings(self, item: dict) -> dict:
Runs anaddb to get the values.
"""
structure = Structure.from_dict(item["abinit_input"]["structure"])
anaddb_input = AnaddbInput.modes_at_qpoint(structure, [0, 0, 0], asr=0, chneut=0)
anaddb_input = AnaddbInput.modes_at_qpoint(
structure, [0, 0, 0], asr=0, chneut=0
)

with tempfile.TemporaryDirectory() as workdir:

Expand Down Expand Up @@ -382,17 +463,25 @@ def get_sum_rule_breakings(self, item: dict) -> dict:
# If the ASR breaking could not be identified. set it to None to signal the
# missing information. This may trigger a warning.
try:
asr_breaking = phbands.asr_breaking(units="cm-1", threshold=0.9, raise_on_no_indices=True)
asr_breaking = phbands.asr_breaking(
units="cm-1", threshold=0.9, raise_on_no_indices=True
)
asr = asr_breaking.absmax_break
except RuntimeError as e:
self.logger.warning("Could not find the ASR breaking for {}. Error: {}".format(item["mp_id"], e))
self.logger.warning(
"Could not find the ASR breaking for {}. Error: {}".format(
item["mp_id"], e
)
)
asr = None

breakings = {"cnsr": cnsr, "asr": asr, "becs_nosymm": becs_val}

return breakings

def run_anaddb(self, ddb_path: str, anaddb_input: AnaddbInput, workdir: str) -> AnaddbTask:
def run_anaddb(
self, ddb_path: str, anaddb_input: AnaddbInput, workdir: str
) -> AnaddbTask:
"""
Runs anaddb. Raise AnaddbError if the calculation couldn't complete
Expand All @@ -404,7 +493,9 @@ def run_anaddb(self, ddb_path: str, anaddb_input: AnaddbInput, workdir: str) ->
An abipy AnaddbTask instance.
"""

task = AnaddbTask.temp_shell_task(anaddb_input, ddb_node=ddb_path, workdir=workdir, manager=self.manager)
task = AnaddbTask.temp_shell_task(
anaddb_input, ddb_node=ddb_path, workdir=workdir, manager=self.manager
)

# Run the task here.
self.logger.debug("Start anaddb for {}".format(ddb_path))
Expand All @@ -420,7 +511,12 @@ def run_anaddb(self, ddb_path: str, anaddb_input: AnaddbInput, workdir: str) ->
return task

def get_properties_anaddb_input(
self, item: dict, bs: bool = True, dos: str = "tetra", lo_to_splitting: bool = True, use_dieflag: bool = True
self,
item: dict,
bs: bool = True,
dos: str = "tetra",
lo_to_splitting: bool = True,
use_dieflag: bool = True,
) -> Tuple[AnaddbInput, Optional[List]]:
"""
creates the AnaddbInput object to calculate the phonon properties.
Expand Down Expand Up @@ -449,7 +545,13 @@ def get_properties_anaddb_input(
inp = AnaddbInput(structure, comment="ANADB input for phonon bands and DOS")

inp.set_vars(
ifcflag=1, ngqpt=np.array(ngqpt), q1shft=q1shft, nqshft=len(q1shft), asr=asr, chneut=chneut, dipdip=dipdip,
ifcflag=1,
ngqpt=np.array(ngqpt),
q1shft=q1shft,
nqshft=len(q1shft),
asr=asr,
chneut=chneut,
dipdip=dipdip,
)

# Parameters for the dos.
Expand All @@ -467,25 +569,35 @@ def get_properties_anaddb_input(
ng2qppa = 500000
dossmear = 1.82e-5 # Ha = 4 cm^-1
ng2qpt = KSampling.automatic_density(structure, kppa=ng2qppa).kpts[0]
inp.set_vars(prtdos=prtdos, dosdeltae=dosdeltae, dossmear=dossmear, ng2qpt=ng2qpt)
inp.set_vars(
prtdos=prtdos, dosdeltae=dosdeltae, dossmear=dossmear, ng2qpt=ng2qpt
)
elif dos is not None:
raise ValueError("Unsupported value of dos.")

# Parameters for the BS
labels_list = None
if bs:
spga = SpacegroupAnalyzer(structure, symprec=self.symprec, angle_tolerance=self.angle_tolerance)
spga = SpacegroupAnalyzer(
structure, symprec=self.symprec, angle_tolerance=self.angle_tolerance
)

spgn = spga.get_space_group_number()
if spgn != item["spacegroup"]["number"]:
raise RuntimeError(
"Parsed specegroup number {} does not match "
"calculation spacegroup {}".format(spgn, item["spacegroup"]["number"])
"calculation spacegroup {}".format(
spgn, item["spacegroup"]["number"]
)
)

hs = HighSymmKpath(structure, symprec=self.symprec, angle_tolerance=self.angle_tolerance)
hs = HighSymmKpath(
structure, symprec=self.symprec, angle_tolerance=self.angle_tolerance
)

qpts, labels_list = hs.get_kpoints(line_density=18, coords_are_cartesian=False)
qpts, labels_list = hs.get_kpoints(
line_density=18, coords_are_cartesian=False
)

n_qpoints = len(qpts)
qph1l = np.zeros((n_qpoints, 4))
Expand Down Expand Up @@ -530,7 +642,9 @@ def get_properties_anaddb_input(
return inp, labels_list

@staticmethod
def get_pmg_bs(phbands: PhononBands, labels_list: List) -> PhononBandStructureSymmLine:
def get_pmg_bs(
phbands: PhononBands, labels_list: List
) -> PhononBandStructureSymmLine:
"""
Generates a PhononBandStructureSymmLine starting from a abipy PhononBands object
Expand Down Expand Up @@ -564,7 +678,9 @@ def get_pmg_bs(phbands: PhononBands, labels_list: List) -> PhononBandStructureSy
displ[i] = phbands._get_non_anal_phdispl(qpts[i + 1])

ph_freqs = np.transpose(ph_freqs) * eV_to_THz
displ = np.transpose(np.reshape(displ, (len(qpts), 3 * n_at, n_at, 3)), (1, 0, 2, 3))
displ = np.transpose(
np.reshape(displ, (len(qpts), 3 * n_at, n_at, 3)), (1, 0, 2, 3)
)

ph_bs_sl = PhononBandStructureSymmLine(
qpoints=qpts,
Expand Down Expand Up @@ -657,7 +773,9 @@ def ensure_indexes(self):
self.phonon_website.ensure_index(self.phonon.key, unique=True)


def get_warnings(asr_break: float, cnsr_break: float, ph_bs: PhononBandStructureSymmLine) -> List[PhononWarnings]:
def get_warnings(
asr_break: float, cnsr_break: float, ph_bs: PhononBandStructureSymmLine
) -> List[PhononWarnings]:
"""
Args:
Expand Down Expand Up @@ -696,7 +814,9 @@ def get_warnings(asr_break: float, cnsr_break: float, ph_bs: PhononBandStructure
return warnings


def get_thermodynamic_properties(ph_dos: CompletePhononDos) -> Tuple[ThermodynamicProperties, VibrationalEnergy]:
def get_thermodynamic_properties(
ph_dos: CompletePhononDos
) -> Tuple[ThermodynamicProperties, VibrationalEnergy]:
"""
Calculates the thermodynamic properties from a phonon DOS
Expand Down

0 comments on commit 83b9b95

Please sign in to comment.