diff --git a/dev_scripts/regen_libxcfunc.py b/dev_scripts/regen_libxcfunc.py index 29e308ccb75..7d114e600b4 100755 --- a/dev_scripts/regen_libxcfunc.py +++ b/dev_scripts/regen_libxcfunc.py @@ -45,28 +45,28 @@ def parse_section(section): return dct -def write_libxc_docs_json(xcfuncs, jpath): +def write_libxc_docs_json(xc_funcs, json_path): """Write json file with libxc metadata to path jpath.""" - xcfuncs = deepcopy(xcfuncs) + xc_funcs = deepcopy(xc_funcs) # Remove XC_FAMILY from Family and XC_ from Kind to make strings more human-readable. - for d in xcfuncs.values(): + for d in xc_funcs.values(): d["Family"] = d["Family"].replace("XC_FAMILY_", "", 1) d["Kind"] = d["Kind"].replace("XC_", "", 1) # Build lightweight version with a subset of keys. - for num, d in xcfuncs.items(): - xcfuncs[num] = {key: d[key] for key in ("Family", "Kind", "References")} + for num, d in xc_funcs.items(): + xc_funcs[num] = {key: d[key] for key in ("Family", "Kind", "References")} # Descriptions are optional for opt in ("Description 1", "Description 2"): desc = d.get(opt) if desc is not None: - xcfuncs[num][opt] = desc + xc_funcs[num][opt] = desc - with open(jpath, mode="w") as file: - json.dump(xcfuncs, file) + with open(json_path, "w") as fh: + json.dump(xc_funcs, fh) - return xcfuncs + return xc_funcs def main(): diff --git a/pymatgen/io/vasp/outputs.py b/pymatgen/io/vasp/outputs.py index 6df0f350cc0..23fceb0cbd4 100644 --- a/pymatgen/io/vasp/outputs.py +++ b/pymatgen/io/vasp/outputs.py @@ -140,6 +140,23 @@ def _vasprun_float(flt: float | str) -> float: raise exc +@dataclass +class KpointOptProps: + """Simple container class to store KPOINTS_OPT data in a separate namespace. Used by Vasprun.""" + + tdos: Dos | None = None + idos: Dos | None = None + pdos: list | None = None + efermi: float | None = None + eigenvalues: dict | None = None + projected_eigenvalues: dict | None = None + projected_magnetisation: np.ndarray | None = None + kpoints: Kpoints | None = None + actual_kpoints: list | None = None + actual_kpoints_weights: list | None = None + dos_has_errors: bool | None = None + + class Vasprun(MSONable): """ Vastly improved cElementTree-based parser for vasprun.xml files. Uses @@ -162,7 +179,7 @@ class Vasprun(MSONable): To access a particular value, you need to do Vasprun.projected_eigenvalues[spin][kpoint index][band index][atom index][orbital_index]. The kpoint, band and atom indices are 0-based (unlike the 1-based indexing in VASP). - projected_magnetisation (numpy.ndarray): Final projected magnetization as a numpy array with the + projected_magnetisation (np.array): Final projected magnetization as a numpy array with the shape (nkpoints, nbands, natoms, norbitals, 3). Where the last axis is the contribution in the 3 Cartesian directions. This attribute is only set if spin-orbit coupling (LSORBIT = True) or non-collinear magnetism (LNONCOLLINEAR = True) is turned on in the INCAR. @@ -174,10 +191,10 @@ class Vasprun(MSONable): real_partyz,real_partxz]],[[imag_partxx,imag_partyy,imag_partzz,imag_partxy, imag_partyz, imag_partxz]]). nionic_steps (int): The total number of ionic steps. This number is always equal to the total number of steps in the actual run even if ionic_step_skip is used. - force_constants (numpy.ndarray): Force constants computed in phonon DFPT run(IBRION = 8). + force_constants (np.array): Force constants computed in phonon DFPT run(IBRION = 8). The data is a 4D numpy array of shape (natoms, natoms, 3, 3). - normalmode_eigenvals (numpy.ndarray): Normal mode frequencies. 1D numpy array of size 3*natoms. - normalmode_eigenvecs (numpy.ndarray): Normal mode eigen vectors. 3D numpy array of shape (3*natoms, natoms, 3). + normalmode_eigenvals (np.array): Normal mode frequencies. 1D numpy array of size 3*natoms. + normalmode_eigenvecs (np.array): Normal mode eigen vectors. 3D numpy array of shape (3*natoms, natoms, 3). md_data (list): Available only for ML MD runs, i.e., INCAR with ML_LMLFF = .TRUE. md_data is a list of dict with the following format: [{'energy': {'e_0_energy': -525.07195568, 'e_fr_energy': -525.07195568, 'e_wo_entrp': -525.07195568, 'kinetic': 3.17809233, 'lattice kinetic': 0.0, 'nosekinetic': 1.323e-5, @@ -192,6 +209,10 @@ class Vasprun(MSONable): 0.04166667, ....]. atomic_symbols (list): List of atomic symbols, e.g., ["Li", "Fe", "Fe", "P", "P", "P"]. potcar_symbols (list): List of POTCAR symbols. e.g., ["PAW_PBE Li 17Jan2003", "PAW_PBE Fe 06Sep2000", ..]. + kpoints_opt_props (object): Object whose attributes are the data from KPOINTS_OPT (if present, + else None). Attributes of the same name have the same format and meaning as Vasprun (or they are + None if absent). Attributes are: tdos, idos, pdos, efermi, eigenvalues, projected_eigenvalues, + projected magnetisation, kpoints, actual_kpoints, actual_kpoints_weights, dos_has_errors. Author: Shyue Ping Ong """ @@ -307,104 +328,162 @@ def _parse(self, stream, parse_dos, parse_eigen, parse_projected_eigen): self.dielectric_data = {} self.other_dielectric = {} self.incar = {} + self.kpoints_opt_props = None + ionic_steps = [] md_data = [] parsed_header = False + in_kpoints_opt = False try: - for _, elem in ET.iterparse(stream): + # When parsing XML, start tags tell us when we have entered a block + # while end tags are when we have actually read the data. + # To know if a particular tag is nested within another block, + # we have to read the start tags (otherwise we will only learn of + # the nesting after we have left the data behind). + # When parsing KPOINTS_OPT data, some of the tags in vasprun.xml + # can only be distinguished from their regular counterparts by + # whether they are nested within another block. This is why we + # must read both start and end tags and have flags to tell us + # when we have entered or left a block. (2024-01-26) + for event, elem in ET.iterparse(stream, events=["start", "end"]): tag = elem.tag - if not parsed_header: - if tag == "generator": - self.generator = self._parse_params(elem) - elif tag == "incar": - self.incar = self._parse_params(elem) - elif tag == "kpoints": - if not hasattr(self, "kpoints"): - self.kpoints, self.actual_kpoints, self.actual_kpoints_weights = self._parse_kpoints(elem) - elif tag == "parameters": - self.parameters = self._parse_params(elem) - elif tag == "structure" and elem.attrib.get("name") == "initialpos": - self.initial_structure = self._parse_structure(elem) - self.final_structure = self.initial_structure - elif tag == "atominfo": - self.atomic_symbols, self.potcar_symbols = self._parse_atominfo(elem) - self.potcar_spec = [ - {"titel": p, "hash": None, "summary_stats": {}} for p in self.potcar_symbols - ] - if tag == "calculation": - parsed_header = True - if not self.parameters.get("LCHIMAG", False): - ionic_steps.append(self._parse_calculation(elem)) - else: - ionic_steps.extend(self._parse_chemical_shielding_calculation(elem)) - elif parse_dos and tag == "dos": - try: - self.tdos, self.idos, self.pdos = self._parse_dos(elem) - self.efermi = self.tdos.efermi - self.dos_has_errors = False - except Exception: - self.dos_has_errors = True - elif parse_eigen and tag == "eigenvalues": - self.eigenvalues = self._parse_eigen(elem) - elif parse_projected_eigen and tag == "projected": - self.projected_eigenvalues, self.projected_magnetisation = self._parse_projected_eigen(elem) - elif tag == "dielectricfunction": - if ( - "comment" not in elem.attrib - or elem.attrib["comment"] - == "INVERSE MACROSCOPIC DIELECTRIC TENSOR (including local field effects in RPA (Hartree))" - ): - if "density" not in self.dielectric_data: - self.dielectric_data["density"] = self._parse_diel(elem) - elif "velocity" not in self.dielectric_data: - # "velocity-velocity" is also named - # "current-current" in OUTCAR - self.dielectric_data["velocity"] = self._parse_diel(elem) + if event == "start": + # The start event tells us when we have entered blocks + if tag == "calculation": + parsed_header = True + elif tag == "eigenvalues_kpoints_opt" or tag == "projected_kpoints_opt": + in_kpoints_opt = True + else: # event == "end": + # The end event happens when we have read a block, so have + # its data. + if not parsed_header: + if tag == "generator": + self.generator = self._parse_params(elem) + elif tag == "incar": + self.incar = self._parse_params(elem) + elif tag == "kpoints": + if not hasattr(self, "kpoints"): + self.kpoints, self.actual_kpoints, self.actual_kpoints_weights = self._parse_kpoints( + elem + ) + elif tag == "parameters": + self.parameters = self._parse_params(elem) + elif tag == "structure" and elem.attrib.get("name") == "initialpos": + self.initial_structure = self._parse_structure(elem) + self.final_structure = self.initial_structure + elif tag == "atominfo": + self.atomic_symbols, self.potcar_symbols = self._parse_atominfo(elem) + self.potcar_spec = [ + {"titel": p, "hash": None, "summary_stats": {}} for p in self.potcar_symbols + ] + if tag == "calculation": + parsed_header = True + if not self.parameters.get("LCHIMAG", False): + ionic_steps.append(self._parse_calculation(elem)) else: - raise NotImplementedError("This vasprun.xml has >2 unlabelled dielectric functions") - else: - comment = elem.attrib["comment"] - # VASP 6+ has labels for the density and current - # derived dielectric constants - if comment == "density-density": - self.dielectric_data["density"] = self._parse_diel(elem) - elif comment == "current-current": - self.dielectric_data["velocity"] = self._parse_diel(elem) + ionic_steps.extend(self._parse_chemical_shielding_calculation(elem)) + elif parse_dos and tag == "dos": + if elem.get("comment") == "kpoints_opt": + if self.kpoints_opt_props is None: + self.kpoints_opt_props = KpointOptProps() + try: + ( + self.kpoints_opt_props.tdos, + self.kpoints_opt_props.idos, + self.kpoints_opt_props.pdos, + ) = self._parse_dos(elem) + self.kpoints_opt_props.efermi = self.kpoints_opt_props.tdos.efermi + self.kpoints_opt_props.dos_has_errors = False + except Exception: + self.kpoints_opt_props.dos_has_errors = True + else: + try: + self.tdos, self.idos, self.pdos = self._parse_dos(elem) + self.efermi = self.tdos.efermi + self.dos_has_errors = False + except Exception: + self.dos_has_errors = True + elif parse_eigen and tag == "eigenvalues" and not in_kpoints_opt: + self.eigenvalues = self._parse_eigen(elem) + elif parse_projected_eigen and tag == "projected" and not in_kpoints_opt: + self.projected_eigenvalues, self.projected_magnetisation = self._parse_projected_eigen(elem) + elif tag == "eigenvalues_kpoints_opt" or tag == "projected_kpoints_opt": + in_kpoints_opt = False + if self.kpoints_opt_props is None: + self.kpoints_opt_props = KpointOptProps() + if parse_eigen: + # projected_kpoints_opt includes occupation information whereas + # eigenvalues_kpoints_opt doesn't. + self.kpoints_opt_props.eigenvalues = self._parse_eigen(elem.find("eigenvalues")) + if tag == "eigenvalues_kpoints_opt": + ( + self.kpoints_opt_props.kpoints, + self.kpoints_opt_props.actual_kpoints, + self.kpoints_opt_props.actual_kpoints_weights, + ) = self._parse_kpoints(elem.find("kpoints")) + elif parse_projected_eigen: # and tag == "projected_kpoints_opt": (implied) + ( + self.kpoints_opt_props.projected_eigenvalues, + self.kpoints_opt_props.projected_magnetisation, + ) = self._parse_projected_eigen(elem) + elif tag == "dielectricfunction": + if ( + "comment" not in elem.attrib + or elem.attrib["comment"] + == "INVERSE MACROSCOPIC DIELECTRIC TENSOR (including local field effects in RPA (Hartree))" + ): + if "density" not in self.dielectric_data: + self.dielectric_data["density"] = self._parse_diel(elem) + elif "velocity" not in self.dielectric_data: + # "velocity-velocity" is also named + # "current-current" in OUTCAR + self.dielectric_data["velocity"] = self._parse_diel(elem) + else: + raise NotImplementedError("This vasprun.xml has >2 unlabelled dielectric functions") else: - self.other_dielectric[comment] = self._parse_diel(elem) - - elif tag == "varray" and elem.attrib.get("name") == "opticaltransitions": - self.optical_transition = np.array(_parse_vasp_array(elem)) - elif tag == "structure" and elem.attrib.get("name") == "finalpos": - self.final_structure = self._parse_structure(elem) - elif tag == "dynmat": - hessian, eigenvalues, eigenvectors = self._parse_dynmat(elem) - # n_atoms is not the total number of atoms, only those for which force constants were calculated - # https://github.com/materialsproject/pymatgen/issues/3084 - n_atoms = len(hessian) // 3 - hessian = np.array(hessian) - self.force_constants = np.zeros((n_atoms, n_atoms, 3, 3), dtype="double") - for ii in range(n_atoms): - for jj in range(n_atoms): - self.force_constants[ii, jj] = hessian[ii * 3 : (ii + 1) * 3, jj * 3 : (jj + 1) * 3] - phonon_eigenvectors = [] - for ev in eigenvectors: - phonon_eigenvectors.append(np.array(ev).reshape(n_atoms, 3)) - self.normalmode_eigenvals = np.array(eigenvalues) - self.normalmode_eigenvecs = np.array(phonon_eigenvectors) - elif self.incar.get("ML_LMLFF"): - if tag == "structure" and elem.attrib.get("name") is None: - md_data.append({}) - md_data[-1]["structure"] = self._parse_structure(elem) - elif tag == "varray" and elem.attrib.get("name") == "forces": - md_data[-1]["forces"] = _parse_vasp_array(elem) - elif tag == "varray" and elem.attrib.get("name") == "stress": - md_data[-1]["stress"] = _parse_vasp_array(elem) - elif tag == "energy": - e_dict = {i.attrib["name"]: float(i.text) for i in elem.findall("i")} - if "kinetic" in e_dict: - md_data[-1]["energy"] = {i.attrib["name"]: float(i.text) for i in elem.findall("i")} + comment = elem.attrib["comment"] + # VASP 6+ has labels for the density and current + # derived dielectric constants + if comment == "density-density": + self.dielectric_data["density"] = self._parse_diel(elem) + elif comment == "current-current": + self.dielectric_data["velocity"] = self._parse_diel(elem) + else: + self.other_dielectric[comment] = self._parse_diel(elem) + + elif tag == "varray" and elem.attrib.get("name") == "opticaltransitions": + self.optical_transition = np.array(_parse_vasp_array(elem)) + elif tag == "structure" and elem.attrib.get("name") == "finalpos": + self.final_structure = self._parse_structure(elem) + elif tag == "dynmat": + hessian, eigenvalues, eigenvectors = self._parse_dynmat(elem) + # n_atoms is not the total number of atoms, only those for which force constants were calculated + # https://github.com/materialsproject/pymatgen/issues/3084 + n_atoms = len(hessian) // 3 + hessian = np.array(hessian) + self.force_constants = np.zeros((n_atoms, n_atoms, 3, 3), dtype="double") + for ii in range(n_atoms): + for jj in range(n_atoms): + self.force_constants[ii, jj] = hessian[ii * 3 : (ii + 1) * 3, jj * 3 : (jj + 1) * 3] + phonon_eigenvectors = [] + for ev in eigenvectors: + phonon_eigenvectors.append(np.array(ev).reshape(n_atoms, 3)) + self.normalmode_eigenvals = np.array(eigenvalues) + self.normalmode_eigenvecs = np.array(phonon_eigenvectors) + elif self.incar.get("ML_LMLFF"): + if tag == "structure" and elem.attrib.get("name") is None: + md_data.append({}) + md_data[-1]["structure"] = self._parse_structure(elem) + elif tag == "varray" and elem.attrib.get("name") == "forces": + md_data[-1]["forces"] = _parse_vasp_array(elem) + elif tag == "varray" and elem.attrib.get("name") == "stress": + md_data[-1]["stress"] = _parse_vasp_array(elem) + elif tag == "energy": + d = {i.attrib["name"]: float(i.text) for i in elem.findall("i")} + if "kinetic" in d: + md_data[-1]["energy"] = {i.attrib["name"]: float(i.text) for i in elem.findall("i")} + except ET.ParseError as exc: if self.exception_on_bad_xml: raise exc @@ -759,6 +838,7 @@ def get_band_structure( efermi: float | Literal["smart"] | None = None, line_mode: bool = False, force_hybrid_mode: bool = False, + ignore_kpoints_opt: bool = False, ) -> BandStructureSymmLine | BandStructure: """Get the band structure as a BandStructure object. @@ -767,7 +847,7 @@ def get_band_structure( the band structure is generated. If none is provided, the code will try to intelligently determine the appropriate KPOINTS file by substituting the - filename of the vasprun.xml with KPOINTS. + filename of the vasprun.xml with KPOINTS (or KPOINTS_OPT). The latter is the default behavior. efermi: The Fermi energy associated with the bandstructure, in eV. By default (None), uses the value reported by VASP in vasprun.xml. To @@ -781,6 +861,8 @@ def get_band_structure( a run along symmetry lines. (Default: False) force_hybrid_mode: Makes it possible to read in self-consistent band structure calculations for every type of functional. (Default: False) + ignore_kpoints_opt: Normally, if KPOINTS_OPT data exists, it has + the band structure data. Set this flag to ignore it. (Default: False) Returns: a BandStructure object (or more specifically a @@ -793,12 +875,20 @@ def get_band_structure( file) (it's not possible to run a non-sc band structure with hybrid functionals). The explicit KPOINTS file needs to have data on the kpoint label as commentary. + + If VASP was run with KPOINTS_OPT, it reads the data from that + file unless told otherwise. This overrides hybrid mode. """ + use_kpoints_opt = not ignore_kpoints_opt and (getattr(self, "kpoints_opt_props", None) is not None) if not kpoints_filename: - kpoints_filename = zpath(os.path.join(os.path.dirname(self.filename), "KPOINTS")) - if kpoints_filename and not os.path.exists(kpoints_filename) and line_mode is True: - raise VaspParseError("KPOINTS not found but needed to obtain band structure along symmetry lines.") - + kpts_path = os.path.join(os.path.dirname(self.filename), "KPOINTS_OPT" if use_kpoints_opt else "KPOINTS") + kpoints_filename = zpath(kpts_path) + if kpoints_filename and not os.path.isfile(kpoints_filename) and line_mode is True: + name = "KPOINTS_OPT" if use_kpoints_opt else "KPOINTS" + raise VaspParseError(f"{name} not found but needed to obtain band structure along symmetry lines.") + + # Note that we're using the Fermi energy of the self-consistent grid + # run even if we're reading bands from KPOINTS_OPT. if efermi == "smart": e_fermi = self.calculate_efermi() elif efermi is None: @@ -811,33 +901,39 @@ def get_band_structure( kpoint_file = Kpoints.from_file(kpoints_filename) lattice_new = Lattice(self.final_structure.lattice.reciprocal_lattice.matrix) - kpoints = [np.array(kpt) for kpt in self.actual_kpoints] + if use_kpoints_opt: + kpoints = [np.array(kpt) for kpt in self.kpoints_opt_props.actual_kpoints] + else: + kpoints = [np.array(kpt) for kpt in self.actual_kpoints] - p_eigenvals: defaultdict[Spin, list] = defaultdict(list) + p_eig_vals: defaultdict[Spin, list] = defaultdict(list) eigenvals: defaultdict[Spin, list] = defaultdict(list) nkpts = len(kpoints) - for spin, v in self.eigenvalues.items(): - v = np.swapaxes(v, 0, 1) - eigenvals[spin] = v[:, :, 0] + if use_kpoints_opt: + eig_vals = self.kpoints_opt_props.eigenvalues + projected_eig_vals = self.kpoints_opt_props.projected_eigenvalues + else: + eig_vals = self.eigenvalues + projected_eig_vals = self.projected_eigenvalues - if self.projected_eigenvalues: - peigen = self.projected_eigenvalues[spin] - # Original axes for self.projected_eigenvalues are kpoints, - # band, ion, orb. + for spin, val in eig_vals.items(): + val = np.swapaxes(val, 0, 1) + eigenvals[spin] = val[:, :, 0] + + if projected_eig_vals: + proj_eig_vals = projected_eig_vals[spin] + # Original axes for self.projected_eigenvalues are kpoints, band, ion, orb. # For BS input, we need band, kpoints, orb, ion. - peigen = np.swapaxes(peigen, 0, 1) # Swap kpoint and band axes - peigen = np.swapaxes(peigen, 2, 3) # Swap ion and orb axes + proj_eig_vals = np.swapaxes(proj_eig_vals, 0, 1) # Swap kpoint and band axes + proj_eig_vals = np.swapaxes(proj_eig_vals, 2, 3) # Swap ion and orb axes - p_eigenvals[spin] = peigen - # for b in range(min_eigenvalues): - # p_eigenvals[spin].append( - # [{Orbital(orb): v for orb, v in enumerate(peigen[b, k])} - # for k in range(nkpts)]) + p_eig_vals[spin] = proj_eig_vals # check if we have an hybrid band structure computation # for this we look at the presence of the LHFCALC tag + # (but hybrid mode is redundant if using kpoints_opt) hybrid_band = False if self.parameters.get("LHFCALC", False) or 0.0 in self.actual_kpoints_weights: hybrid_band = True @@ -847,7 +943,7 @@ def get_band_structure( if line_mode: labels_dict = {} - if hybrid_band or force_hybrid_mode: + if (hybrid_band or force_hybrid_mode) and not use_kpoints_opt: start_bs_index = 0 for i in range(len(self.actual_kpoints)): if self.actual_kpoints_weights[i] == 0.0: @@ -862,15 +958,13 @@ def get_band_structure( kpoints = kpoints[start_bs_index:nkpts] up_eigen = [eigenvals[Spin.up][i][start_bs_index:nkpts] for i in range(nbands)] if self.projected_eigenvalues: - p_eigenvals[Spin.up] = [p_eigenvals[Spin.up][i][start_bs_index:nkpts] for i in range(nbands)] + p_eig_vals[Spin.up] = [p_eig_vals[Spin.up][i][start_bs_index:nkpts] for i in range(nbands)] if self.is_spin: down_eigen = [eigenvals[Spin.down][i][start_bs_index:nkpts] for i in range(nbands)] eigenvals[Spin.up] = up_eigen eigenvals[Spin.down] = down_eigen if self.projected_eigenvalues: - p_eigenvals[Spin.down] = [ - p_eigenvals[Spin.down][i][start_bs_index:nkpts] for i in range(nbands) - ] + p_eig_vals[Spin.down] = [p_eig_vals[Spin.down][i][start_bs_index:nkpts] for i in range(nbands)] else: eigenvals[Spin.up] = up_eigen else: @@ -889,7 +983,7 @@ def get_band_structure( e_fermi, labels_dict, structure=self.final_structure, - projections=p_eigenvals, + projections=p_eig_vals, ) return BandStructure( kpoints, # type: ignore[arg-type] @@ -897,7 +991,7 @@ def get_band_structure( lattice_new, e_fermi, structure=self.final_structure, - projections=p_eigenvals, # type: ignore[arg-type] + projections=p_eig_vals, # type: ignore[arg-type] ) @property @@ -1117,13 +1211,24 @@ def as_dict(self): } actual_kpts = [ { - "abc": list(self.actual_kpoints[i]), - "weight": self.actual_kpoints_weights[i], + "abc": list(self.actual_kpoints[idx]), + "weight": self.actual_kpoints_weights[idx], } - for i in range(len(self.actual_kpoints)) + for idx in range(len(self.actual_kpoints)) ] vin["kpoints"]["actual_points"] = actual_kpts vin["nkpoints"] = len(actual_kpts) + if kpt_opt_props := getattr(self, "kpoints_opt_props", None): + vin["kpoints_opt"] = kpt_opt_props.kpoints.as_dict() + actual_kpts = [ + { + "abc": list(kpt_opt_props.actual_kpoints[idx]), + "weight": kpt_opt_props.actual_kpoints_weights[idx], + } + for idx in range(len(kpt_opt_props.actual_kpoints)) + ] + vin["kpoints_opt"]["actual_kpoints"] = actual_kpts + vin["nkpoints_opt"] = len(actual_kpts) vin["potcar"] = [s.split(" ")[1] for s in self.potcar_symbols] vin["potcar_spec"] = self.potcar_spec vin["potcar_type"] = [s.split(" ")[0] for s in self.potcar_symbols] @@ -1164,6 +1269,21 @@ def as_dict(self): if self.projected_magnetisation is not None: vout["projected_magnetisation"] = self.projected_magnetisation.tolist() + if kpt_opt_props and kpt_opt_props.eigenvalues: + eigen = {str(spin): v.tolist() for spin, v in kpt_opt_props.eigenvalues.items()} + vout["eigenvalues_kpoints_opt"] = eigen + # TODO implement kpoints_opt eigenvalue_band_proprties. + # (gap, cbm, vbm, is_direct) = self.eigenvalue_band_properties + # vout.update({"bandgap": gap, "cbm": cbm, "vbm": vbm, "is_gap_direct": is_direct}) + + if kpt_opt_props.projected_eigenvalues: + vout["projected_eigenvalues_kpoints_opt"] = { + str(spin): v.tolist() for spin, v in kpt_opt_props.projected_eigenvalues.items() + } + + if kpt_opt_props.projected_magnetisation is not None: + vout["projected_magnetisation_kpoints_opt"] = kpt_opt_props.projected_magnetisation.tolist() + vout["epsilon_static"] = self.epsilon_static vout["epsilon_static_wolfe"] = self.epsilon_static_wolfe vout["epsilon_ionic"] = self.epsilon_ionic @@ -1477,36 +1597,71 @@ def __init__( with zopen(filename, mode="rt") as file: self.efermi = None parsed_header = False + in_kpoints_opt = False self.eigenvalues = self.projected_eigenvalues = None - for _, elem in ET.iterparse(file): + self.kpoints_opt_props = None + for event, elem in ET.iterparse(file, events=["start", "end"]): tag = elem.tag - if not parsed_header: - if tag == "generator": - self.generator = self._parse_params(elem) - elif tag == "incar": - self.incar = self._parse_params(elem) - elif tag == "kpoints": - ( - self.kpoints, - self.actual_kpoints, - self.actual_kpoints_weights, - ) = self._parse_kpoints(elem) - elif tag == "parameters": - self.parameters = self._parse_params(elem) - elif tag == "atominfo": - self.atomic_symbols, self.potcar_symbols = self._parse_atominfo(elem) - self.potcar_spec = [ - {"titel": p, "hash": None, "summary_stats": {}} for p in self.potcar_symbols - ] - parsed_header = True - elif tag == "i" and elem.attrib.get("name") == "efermi": - self.efermi = float(elem.text) - elif tag == "eigenvalues": - self.eigenvalues = self._parse_eigen(elem) - elif parse_projected_eigen and tag == "projected": - self.projected_eigenvalues, self.projected_magnetisation = self._parse_projected_eigen(elem) - elif tag == "structure" and elem.attrib.get("name") == "finalpos": - self.final_structure = self._parse_structure(elem) + if event == "start": + # The start event tells us when we have entered blocks + if ( + tag == "eigenvalues_kpoints_opt" + or tag == "projected_kpoints_opt" + or (tag == "dos" and elem.attrib.get("comment") == "kpoints_opt") + ): + in_kpoints_opt = True + else: # if event == "end": + if not parsed_header: + if tag == "generator": + self.generator = self._parse_params(elem) + elif tag == "incar": + self.incar = self._parse_params(elem) + elif tag == "kpoints": + ( + self.kpoints, + self.actual_kpoints, + self.actual_kpoints_weights, + ) = self._parse_kpoints(elem) + elif tag == "parameters": + self.parameters = self._parse_params(elem) + elif tag == "atominfo": + self.atomic_symbols, self.potcar_symbols = self._parse_atominfo(elem) + self.potcar_spec = [ + {"titel": p, "hash": None, "summary_stats": {}} for p in self.potcar_symbols + ] + parsed_header = True + elif tag == "i" and elem.attrib.get("name") == "efermi": + if in_kpoints_opt: + if self.kpoints_opt_props is None: + self.kpoints_opt_props = KpointOptProps() + self.kpoints_opt_props.efermi = float(elem.text) + in_kpoints_opt = False + else: + self.efermi = float(elem.text) + elif tag == "eigenvalues" and not in_kpoints_opt: + self.eigenvalues = self._parse_eigen(elem) + elif parse_projected_eigen and tag == "projected" and not in_kpoints_opt: + self.projected_eigenvalues, self.projected_magnetisation = self._parse_projected_eigen(elem) + elif tag == "eigenvalues_kpoints_opt" or tag == "projected_kpoints_opt": + if self.kpoints_opt_props is None: + self.kpoints_opt_props = KpointOptProps() + in_kpoints_opt = False + # projected_kpoints_opt includes occupation information whereas + # eigenvalues_kpoints_opt doesn't. + self.kpoints_opt_props.eigenvalues = self._parse_eigen(elem.find("eigenvalues")) + if tag == "eigenvalues_kpoints_opt": + ( + self.kpoints_opt_props.kpoints, + self.kpoints_opt_props.actual_kpoints, + self.kpoints_opt_props.actual_kpoints_weights, + ) = self._parse_kpoints(elem.find("kpoints")) + elif parse_projected_eigen: # and tag == "projected_kpoints_opt": (implied) + ( + self.kpoints_opt_props.projected_eigenvalues, + self.kpoints_opt_props.projected_magnetisation, + ) = self._parse_projected_eigen(elem) + elif tag == "structure" and elem.attrib.get("name") == "finalpos": + self.final_structure = self._parse_structure(elem) self.vasp_version = self.generator["version"] if parse_potcar_file: self.update_potcar_spec(parse_potcar_file) @@ -1544,6 +1699,17 @@ def as_dict(self): for i in range(len(self.actual_kpoints)) ] vin["kpoints"]["actual_points"] = actual_kpts + if kpt_opt_props := getattr(self, "kpoints_opt_props", None): + vin["kpoints_opt"] = kpt_opt_props.kpoints.as_dict() + actual_kpts = [ + { + "abc": list(kpt_opt_props.actual_kpoints[idx]), + "weight": kpt_opt_props.actual_kpoints_weights[idx], + } + for idx in range(len(kpt_opt_props.actual_kpoints)) + ] + vin["kpoints_opt"]["actual_kpoints"] = actual_kpts + vin["nkpoints_opt"] = len(actual_kpts) vin["potcar"] = [s.split(" ")[1] for s in self.potcar_symbols] vin["potcar_spec"] = self.potcar_spec vin["potcar_type"] = [s.split(" ")[0] for s in self.potcar_symbols] @@ -1570,6 +1736,18 @@ def as_dict(self): peigen[kpoint_index][str(spin)] = vv vout["projected_eigenvalues"] = peigen + if kpt_opt_props and kpt_opt_props.eigenvalues: + eigen = {str(spin): v.tolist() for spin, v in kpt_opt_props.eigenvalues.items()} + vout["eigenvalues_kpoints_opt"] = eigen + # TODO implement kpoints_opt eigenvalue_band_proprties. + # (gap, cbm, vbm, is_direct) = self.eigenvalue_band_properties + # vout.update({"bandgap": gap, "cbm": cbm, "vbm": vbm, "is_gap_direct": is_direct}) + + if kpt_opt_props.projected_eigenvalues: + vout["projected_eigenvalues_kpoints_opt"] = { + str(spin): v.tolist() for spin, v in kpt_opt_props.projected_eigenvalues.items() + } + dct["output"] = vout return jsanitize(dct, strict=True) @@ -1591,7 +1769,7 @@ class Outcar: chemical_shielding (dict): Chemical shielding on each ion as a dictionary with core and valence contributions. unsym_cs_tensor (list): Unsymmetrized chemical shielding tensor matrixes on each ion as a list. e.g., [[[sigma11, sigma12, sigma13], [sigma21, sigma22, sigma23], [sigma31, sigma32, sigma33]], ...] - cs_g0_contribution (numpy.ndarray): G=0 contribution to chemical shielding. 2D rank 3 matrix. + cs_g0_contribution (np.array): G=0 contribution to chemical shielding. 2D rank 3 matrix. cs_core_contribution (dict): Core contribution to chemical shielding. dict. e.g., {'Mg': -412.8, 'C': -200.5, 'O': -271.1} efg (tuple): Electric Field Gradient (EFG) tensor on each ion as a tuple of dict, e.g., @@ -1602,12 +1780,12 @@ class Outcar: is_stopped (bool): True if OUTCAR is from a stopped run (using STOPCAR, see VASP Manual). run_stats (dict): Various useful run stats as a dict including "System time (sec)", "Total CPU time used (sec)", "Elapsed time (sec)", "Maximum memory used (kb)", "Average memory used (kb)", "User time (sec)", "cores". - elastic_tensor (numpy.ndarray): Total elastic moduli (Kbar) is given in a 6x6 array matrix. - drift (numpy.ndarray): Total drift for each step in eV/Atom. + elastic_tensor (np.array): Total elastic moduli (Kbar) is given in a 6x6 array matrix. + drift (np.array): Total drift for each step in eV/Atom. ngf (tuple): Dimensions for the Augmentation grid. - sampling_radii (numpy.ndarray): Size of the sampling radii in VASP for the test charges for the electrostatic + sampling_radii (np.array): Size of the sampling radii in VASP for the test charges for the electrostatic potential at each atom. Total array size is the number of elements present in the calculation. - electrostatic_potential (numpy.ndarray): Average electrostatic potential at each atomic position in order of + electrostatic_potential (np.array): Average electrostatic potential at each atomic position in order of the atoms in POSCAR. final_energy_contribs (dict): Individual contributions to the total final energy as a dictionary. Include contributions from keys, e.g.: @@ -3564,7 +3742,7 @@ class Procar: data (dict): The PROCAR data of the form below. It should VASP uses 1-based indexing, but all indices are converted to 0-based here. { spin: nd.array accessed with (k-point index, band index, ion index, orbital index) } - weights (numpy.ndarray): The weights associated with each k-point as an nd.array of length nkpoints. + weights (np.array): The weights associated with each k-point as an nd.array of length nkpoints. phase_factors (dict): Phase factors, where present (e.g. LORBIT = 12). A dict of the form: { spin: complex nd.array accessed with (k-point index, band index, ion index, orbital index) } nbands (int): Number of bands. @@ -4194,10 +4372,10 @@ class Wavecar: nb (int): Number of bands per k-point. encut (float): Energy cutoff (used to define G_{cut}). efermi (float): Fermi energy. - a (numpy.ndarray): Primitive lattice vectors of the cell (e.g. a_1 = self.a[0, :]). - b (numpy.ndarray): Reciprocal lattice vectors of the cell (e.g. b_1 = self.b[0, :]). + a (np.array): Primitive lattice vectors of the cell (e.g. a_1 = self.a[0, :]). + b (np.array): Reciprocal lattice vectors of the cell (e.g. b_1 = self.b[0, :]). vol (float): The volume of the unit cell in real space. - kpoints (numpy.ndarray): The list of k-points read from the WAVECAR file. + kpoints (np.array): The list of k-points read from the WAVECAR file. band_energy (list): The list of band eigenenergies (and corresponding occupancies) for each kpoint, where the first index corresponds to the index of the k-point (e.g. self.band_energy[kp]). Gpoints (list): The list of generated G-points for each k-point (a double list), which diff --git a/tests/apps/borg/test_queen.py b/tests/apps/borg/test_queen.py index 4eea03b7bb7..c6be59aa7be 100644 --- a/tests/apps/borg/test_queen.py +++ b/tests/apps/borg/test_queen.py @@ -11,9 +11,6 @@ __author__ = "Shyue Ping Ong" __copyright__ = "Copyright 2012, The Materials Project" -__version__ = "0.1" -__maintainer__ = "Shyue Ping Ong" -__email__ = "shyue@mit.edu" __date__ = "Mar 18, 2012" @@ -22,7 +19,7 @@ def test_get_data(self): drone = VaspToComputedEntryDrone() self.queen = BorgQueen(drone, TEST_FILES_DIR, 1) data = self.queen.get_data() - assert len(data) == 15 + assert len(data) == 16 def test_load_data(self): drone = VaspToComputedEntryDrone() diff --git a/tests/files/kpoints_opt/KPOINTS b/tests/files/kpoints_opt/KPOINTS new file mode 100644 index 00000000000..b381769defe --- /dev/null +++ b/tests/files/kpoints_opt/KPOINTS @@ -0,0 +1,4 @@ +Fully automatic kpoint scheme +0 +Automatic +25 diff --git a/tests/files/kpoints_opt/KPOINTS_OPT b/tests/files/kpoints_opt/KPOINTS_OPT new file mode 100644 index 00000000000..fb17416f794 --- /dev/null +++ b/tests/files/kpoints_opt/KPOINTS_OPT @@ -0,0 +1,33 @@ +Line_mode KPOINTS file +10 +Line_mode +Reciprocal +0.0 0.0 0.0 ! \Gamma +0.5 0.0 0.5 ! X + +0.5 0.0 0.5 ! X +0.5 0.25 0.75 ! W + +0.5 0.25 0.75 ! W +0.375 0.375 0.75 ! K + +0.375 0.375 0.75 ! K +0.0 0.0 0.0 ! \Gamma + +0.0 0.0 0.0 ! \Gamma +0.5 0.5 0.5 ! L + +0.5 0.5 0.5 ! L +0.625 0.25 0.625 ! U + +0.625 0.25 0.625 ! U +0.5 0.25 0.75 ! W + +0.5 0.25 0.75 ! W +0.5 0.5 0.5 ! L + +0.5 0.5 0.5 ! L +0.375 0.375 0.75 ! K + +0.625 0.25 0.625 ! U +0.5 0.0 0.5 ! X diff --git a/tests/files/kpoints_opt/vasprun.xml.gz b/tests/files/kpoints_opt/vasprun.xml.gz new file mode 100644 index 00000000000..5a578a1c16b Binary files /dev/null and b/tests/files/kpoints_opt/vasprun.xml.gz differ diff --git a/tests/io/vasp/test_outputs.py b/tests/io/vasp/test_outputs.py index 890b9604606..7c9fe266be8 100644 --- a/tests/io/vasp/test_outputs.py +++ b/tests/io/vasp/test_outputs.py @@ -18,6 +18,7 @@ from pymatgen.core import Element from pymatgen.core.lattice import Lattice from pymatgen.core.structure import Structure +from pymatgen.electronic_structure.bandstructure import BandStructureSymmLine from pymatgen.electronic_structure.core import Magmom, Orbital, OrbitalType, Spin from pymatgen.entries.compatibility import MaterialsProjectCompatibility from pymatgen.io.vasp.inputs import Kpoints, Poscar, Potcar @@ -47,6 +48,8 @@ except ImportError: h5py = None +kpts_opt_vrun_path = f"{TEST_FILES_DIR}/kpoints_opt/vasprun.xml.gz" + class TestVasprun(PymatgenTest): def test_vasprun_ml(self): @@ -72,7 +75,9 @@ def test_vasprun_md(self): assert vasp_run.converged_ionic def test_bad_random_seed(self): - _ = Vasprun(f"{TEST_FILES_DIR}/vasprun.bad_random_seed.xml.gz") + vasp_run = Vasprun(f"{TEST_FILES_DIR}/vasprun.bad_random_seed.xml.gz") + assert vasp_run.incar["ISMEAR"] == 0 + assert vasp_run.incar["RANDOM_SEED"] is None def test_multiple_dielectric(self): vasp_run = Vasprun(f"{TEST_FILES_DIR}/vasprun.GW0.xml.gz") @@ -704,6 +709,59 @@ def test_eigenvalue_band_properties_separate_spins(self): assert props[3][0] assert props[3][1] + def test_kpoints_opt(self): + vasp_run = Vasprun(kpts_opt_vrun_path, parse_projected_eigen=True) + # This calculation was run using KPOINTS_OPT + kpt_opt_props = vasp_run.kpoints_opt_props + # Check the k-points were read correctly. + assert len(vasp_run.actual_kpoints) == 10 + assert len(vasp_run.actual_kpoints_weights) == 10 + assert len(kpt_opt_props.actual_kpoints) == 100 + assert len(kpt_opt_props.actual_kpoints_weights) == 100 + # Check the eigenvalues were read correctly. + assert vasp_run.eigenvalues[Spin.up].shape == (10, 24, 2) + assert kpt_opt_props.eigenvalues[Spin.up].shape == (100, 24, 2) + assert vasp_run.eigenvalues[Spin.up][0, 0, 0] == approx(-6.1471) + assert kpt_opt_props.eigenvalues[Spin.up][0, 0, 0] == approx(-6.1536) + # Check the projected eigenvalues were read correctly + assert vasp_run.projected_eigenvalues[Spin.up].shape == (10, 24, 8, 9) + assert kpt_opt_props.projected_eigenvalues[Spin.up].shape == (100, 24, 8, 9) + assert vasp_run.projected_eigenvalues[Spin.up][0, 1, 0, 0] == approx(0.0492) + # I think these zeroes are a bug in VASP (maybe my VASP) transcribing from PROCAR_OPT to vasprun.xml + # No matter. The point of the parser is to read what's in the file. + assert kpt_opt_props.projected_eigenvalues[Spin.up][0, 1, 0, 0] == approx(0.0000) + # Test as_dict + vasp_run_dct = vasp_run.as_dict() + assert vasp_run_dct["input"]["nkpoints_opt"] == 100 + assert vasp_run_dct["input"]["nkpoints"] == 10 + assert vasp_run_dct["output"]["eigenvalues_kpoints_opt"]["1"][0][0][0] == approx(-6.1536) + + def test_kpoints_opt_band_structure(self): + vasp_run = Vasprun(kpts_opt_vrun_path, parse_potcar_file=False, parse_projected_eigen=True) + bs = vasp_run.get_band_structure(f"{TEST_FILES_DIR}/kpoints_opt/KPOINTS_OPT") + assert isinstance(bs, BandStructureSymmLine) + cbm = bs.get_cbm() + vbm = bs.get_vbm() + assert cbm["kpoint_index"] == [38], "wrong cbm kpoint index" + assert cbm["energy"] == approx(6.4394), "wrong cbm energy" + assert cbm["band_index"] == {Spin.up: [16], Spin.down: [16]}, "wrong cbm bands" + assert vbm["kpoint_index"] == [0, 39, 40] + assert vbm["energy"] == approx(5.7562), "wrong vbm energy" + assert vbm["band_index"] == {Spin.down: [13, 14, 15], Spin.up: [13, 14, 15]}, "wrong vbm bands" + vbm_kp_label = vbm["kpoint"].label + assert vbm["kpoint"].label == "\\Gamma", f"Unpexpected {vbm_kp_label=}" + cmb_kp_label = cbm["kpoint"].label + assert cmb_kp_label is None, f"Unpexpected {cmb_kp_label=}" + # Test projection + projected = bs.get_projection_on_elements() + assert np.isnan(projected[Spin.up][0][0]["Si"]) + # Due to some error in my VASP, the transcription of PROCAR_OPT into + # vasprun.xml is filled to the brim with errors in the projections. + # At some point we might get a healthier vasprun.xml, but the point here + # is to test the parser, not VASP. + projected = bs.get_projections_on_elements_and_orbitals({"Si": ["s"]}) + assert projected[Spin.up][0][58]["Si"]["s"] == -0.0271 + def test_parse_potcar_cwd_relative(self): # Test to ensure that common use cases of vasprun parsing work # in the current working directory using relative paths, @@ -1322,6 +1380,35 @@ def test_get_band_structure(self): dct = vasprun.as_dict() assert "eigenvalues" in dct["output"] + def test_kpoints_opt(self): + vasp_run = BSVasprun(kpts_opt_vrun_path, parse_potcar_file=False, parse_projected_eigen=True) + bs = vasp_run.get_band_structure(f"{TEST_FILES_DIR}/kpoints_opt/KPOINTS_OPT") + assert isinstance(bs, BandStructureSymmLine) + cbm = bs.get_cbm() + vbm = bs.get_vbm() + assert cbm["kpoint_index"] == [38], "wrong cbm kpoint index" + assert cbm["energy"] == approx(6.4394), "wrong cbm energy" + assert cbm["band_index"] == {Spin.up: [16], Spin.down: [16]}, "wrong cbm bands" + # Strangely, when I call with parse_projected_eigen, it gives empty Spin.down, + # but without parse_projected_eigen it does not give it. + # So at one point it called the empty key. + assert vbm["kpoint_index"] == [0, 39, 40] + assert vbm["energy"] == approx(5.7562), "wrong vbm energy" + assert vbm["band_index"] == {Spin.down: [13, 14, 15], Spin.up: [13, 14, 15]}, "wrong vbm bands" + assert vbm["kpoint"].label == "\\Gamma", "wrong vbm label" + assert cbm["kpoint"].label is None, "wrong cbm label" + # Test projection + projected = bs.get_projection_on_elements() + assert np.isnan(projected[Spin.up][0][0]["Si"]) + # Due to some error in my VASP, the transcription of PROCAR_OPT into + # vasprun.xml is filled to the brim with errors in the projections. + # At some point we might get a healthier vasprun.xml, but the point here + # is to test the parser, not VASP. + projected = bs.get_projections_on_elements_and_orbitals({"Si": ["s"]}) + assert projected[Spin.up][0][58]["Si"]["s"] == -0.0271 + vrun_dct = vasp_run.as_dict() + assert {*vrun_dct["output"]} >= {"eigenvalues", "eigenvalues_kpoints_opt"} + class TestOszicar(PymatgenTest): def test_init(self):