Source code for pymatgen.io.pwscf

# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.


"""
This module implements input and output processing from PWSCF.
"""


import re

from monty.io import zopen

from monty.re import regrep
from collections import defaultdict

from pymatgen.core.periodic_table import Element
from pymatgen.core.lattice import Lattice
from pymatgen.core.structure import Structure
from pymatgen.util.io_utils import clean_lines


__author__ = "Shyue Ping Ong"
__copyright__ = "Copyright 2012, The Materials Virtual Lab"
__version__ = "0.1"
__maintainer__ = "Shyue Ping Ong"
__email__ = "ongsp@ucsd.edu"
__date__ = "3/27/15"


[docs]class PWInput: """ Base input file class. Right now, only supports no symmetry and is very basic. """ def __init__(self, structure, pseudo=None, control=None, system=None, electrons=None, ions=None, cell=None, kpoints_mode="automatic", kpoints_grid=(1, 1, 1), kpoints_shift=(0, 0, 0)): """ Initializes a PWSCF input file. Args: structure (Structure): Input structure. For spin-polarized calculation, properties (e.g. {"starting_magnetization": -0.5, "pseudo": "Mn.pbe-sp-van.UPF"}) on each site is needed instead of pseudo (dict). pseudo (dict): A dict of the pseudopotentials to use. Default to None. control (dict): Control parameters. Refer to official PWSCF doc on supported parameters. Default to {"calculation": "scf"} system (dict): System parameters. Refer to official PWSCF doc on supported parameters. Default to None, which means {}. electrons (dict): Electron parameters. Refer to official PWSCF doc on supported parameters. Default to None, which means {}. ions (dict): Ions parameters. Refer to official PWSCF doc on supported parameters. Default to None, which means {}. cell (dict): Cell parameters. Refer to official PWSCF doc on supported parameters. Default to None, which means {}. kpoints_mode (str): Kpoints generation mode. Default to automatic. kpoints_grid (sequence): The kpoint grid. Default to (1, 1, 1). kpoints_shift (sequence): The shift for the kpoints. Defaults to (0, 0, 0). """ self.structure = structure sections = {} sections["control"] = control or {"calculation": "scf"} sections["system"] = system or {} sections["electrons"] = electrons or {} sections["ions"] = ions or {} sections["cell"] = cell or {} if pseudo is None: for site in structure: try: site.properties['pseudo'] except KeyError: raise PWInputError("Missing %s in pseudo specification!" % site) else: for species in self.structure.composition.keys(): if species.symbol not in pseudo: raise PWInputError("Missing %s in pseudo specification!" % species.symbol) self.pseudo = pseudo self.sections = sections self.kpoints_mode = kpoints_mode self.kpoints_grid = kpoints_grid self.kpoints_shift = kpoints_shift def __str__(self): out = [] site_descriptions = {} if self.pseudo is not None: site_descriptions = self.pseudo else: c = 1 for site in self.structure: name = None for k, v in site_descriptions.items(): if site.properties == v: name = k if name is None: name = site.specie.symbol + str(c) site_descriptions[name] = site.properties c += 1 def to_str(v): if isinstance(v, str): return "'%s'" % v elif isinstance(v, float): return "%s" % str(v).replace("e", "d") elif isinstance(v, bool): if v: return ".TRUE." else: return ".FALSE." return v for k1 in ["control", "system", "electrons", "ions", "cell"]: v1 = self.sections[k1] out.append("&%s" % k1.upper()) sub = [] for k2 in sorted(v1.keys()): if isinstance(v1[k2], list): n = 1 for l in v1[k2][:len(site_descriptions)]: sub.append(" %s(%d) = %s" % (k2, n, to_str(v1[k2][n - 1]))) n += 1 else: sub.append(" %s = %s" % (k2, to_str(v1[k2]))) if k1 == "system": if 'ibrav' not in self.sections[k1]: sub.append(" ibrav = 0") if 'nat' not in self.sections[k1]: sub.append(" nat = %d" % len(self.structure)) if 'ntyp' not in self.sections[k1]: sub.append(" ntyp = %d" % len(site_descriptions)) sub.append("/") out.append(",\n".join(sub)) out.append("ATOMIC_SPECIES") for k, v in sorted(site_descriptions.items(), key=lambda i: i[0]): e = re.match(r"[A-Z][a-z]?", k).group(0) if self.pseudo is not None: p = v else: p = v['pseudo'] out.append(" %s %.4f %s" % (k, Element(e).atomic_mass, p)) out.append("ATOMIC_POSITIONS crystal") if self.pseudo is not None: for site in self.structure: out.append(" %s %.6f %.6f %.6f" % (site.specie.symbol, site.a, site.b, site.c)) else: for site in self.structure: name = None for k, v in sorted(site_descriptions.items(), key=lambda i: i[0]): if v == site.properties: name = k out.append(" %s %.6f %.6f %.6f" % (name, site.a, site.b, site.c)) out.append("K_POINTS %s" % self.kpoints_mode) kpt_str = ["%s" % i for i in self.kpoints_grid] kpt_str.extend(["%s" % i for i in self.kpoints_shift]) out.append(" %s" % " ".join(kpt_str)) out.append("CELL_PARAMETERS angstrom") for vec in self.structure.lattice.matrix: out.append(" %f %f %f" % (vec[0], vec[1], vec[2])) return "\n".join(out)
[docs] def as_dict(self): """ Create a dictionary representation of a PWInput object Returns: dict """ pwinput_dict = {'structure': self.structure.as_dict(), 'pseudo': self.pseudo, 'sections': self.sections, 'kpoints_mode': self.kpoints_mode, 'kpoints_grid': self.kpoints_grid, 'kpoints_shift': self.kpoints_shift} return pwinput_dict
[docs] @classmethod def from_dict(cls, pwinput_dict): """ Load a PWInput object from a dictionary. Args: pwinput_dict (dict): dictionary with PWInput data Returns: PWInput object """ pwinput = cls(structure=Structure.from_dict(pwinput_dict['structure']), pseudo=pwinput_dict['pseudo'], control=pwinput_dict['sections']['control'], system=pwinput_dict['sections']['system'], electrons=pwinput_dict['sections']['electrons'], ions=pwinput_dict['sections']['ions'], cell=pwinput_dict['sections']['cell'], kpoints_mode=pwinput_dict['kpoints_mode'], kpoints_grid=pwinput_dict['kpoints_grid'], kpoints_shift=pwinput_dict['kpoints_shift']) return pwinput
[docs] def write_file(self, filename): """ Write the PWSCF input file. Args: filename (str): The string filename to output to. """ with open(filename, "w") as f: f.write(self.__str__())
[docs] @staticmethod def from_file(filename): """ Reads an PWInput object from a file. Args: filename (str): Filename for file Returns: PWInput object """ with zopen(filename, "rt") as f: return PWInput.from_string(f.read())
[docs] @staticmethod def from_string(string): """ Reads an PWInput object from a string. Args: string (str): PWInput string Returns: PWInput object """ lines = list(clean_lines(string.splitlines())) def input_mode(line): if line[0] == "&": return ("sections", line[1:].lower()) elif "ATOMIC_SPECIES" in line: return ("pseudo",) elif "K_POINTS" in line: return ("kpoints", line.split("{")[1][:-1]) elif "CELL_PARAMETERS" in line or "ATOMIC_POSITIONS" in line: return ("structure", line.split("{")[1][:-1]) elif line == "/": return None else: return mode sections = {"control": {}, "system": {}, "electrons": {}, "ions": {}, "cell": {}} pseudo = {} pseudo_index = 0 lattice = [] species = [] coords = [] structure = None site_properties = {"pseudo": []} mode = None for line in lines: mode = input_mode(line) if mode is None: pass elif mode[0] == "sections": section = mode[1] m = re.match(r'(\w+)\(?(\d*?)\)?\s*=\s*(.*)', line) if m: key = m.group(1).strip() key_ = m.group(2).strip() val = m.group(3).strip() if key_ != "": if sections[section].get(key, None) is None: val_ = [0.0] * 20 # MAX NTYP DEFINITION val_[int(key_) - 1] = PWInput.proc_val(key, val) sections[section][key] = val_ site_properties[key] = [] else: sections[section][key][int(key_) - 1] = PWInput.proc_val(key, val) else: sections[section][key] = PWInput.proc_val(key, val) elif mode[0] == "pseudo": m = re.match(r'(\w+)\s+(\d*.\d*)\s+(.*)', line) if m: pseudo[m.group(1).strip()] = {} pseudo[m.group(1).strip()]["index"] = pseudo_index pseudo[m.group(1).strip()]["pseudopot"] = m.group(3).strip() pseudo_index += 1 elif mode[0] == "kpoints": m = re.match(r'(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)', line) if m: kpoints_grid = (int(m.group(1)), int(m.group(2)), int(m.group(3))) kpoints_shift = (int(m.group(4)), int(m.group(5)), int(m.group(6))) else: kpoints_mode = mode[1] elif mode[0] == "structure": m_l = re.match(r'(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)', line) m_p = re.match(r'(\w+)\s+(-?\d+\.\d*)\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)', line) if m_l: lattice += [float(m_l.group(1)), float(m_l.group(2)), float(m_l.group(3))] elif m_p: site_properties["pseudo"].append(pseudo[m_p.group(1)]["pseudopot"]) species += [pseudo[m_p.group(1)]["pseudopot"].split(".")[0]] coords += [[float(m_p.group(2)), float(m_p.group(3)), float(m_p.group(4))]] for k, v in site_properties.items(): if k != "pseudo": site_properties[k].append(sections['system'][k][pseudo[m_p.group(1)]["index"]]) if mode[1] == "angstrom": coords_are_cartesian = True elif mode[1] == "crystal": coords_are_cartesian = False structure = Structure(Lattice(lattice), species, coords, coords_are_cartesian=coords_are_cartesian, site_properties=site_properties) return PWInput(structure=structure, control=sections["control"], system=sections["system"], electrons=sections["electrons"], ions=sections["ions"], cell=sections["cell"], kpoints_mode=kpoints_mode, kpoints_grid=kpoints_grid, kpoints_shift=kpoints_shift)
[docs] def proc_val(key, val): """ Static helper method to convert PWINPUT parameters to proper type, e.g., integers, floats, etc. Args: key: PWINPUT parameter key val: Actual value of PWINPUT parameter. """ float_keys = ('etot_conv_thr', 'forc_conv_thr', 'conv_thr', 'Hubbard_U', 'Hubbard_J0', 'defauss', 'starting_magnetization',) int_keys = ('nstep', 'iprint', 'nberrycyc', 'gdir', 'nppstr', 'ibrav', 'nat', 'ntyp', 'nbnd', 'nr1', 'nr2', 'nr3', 'nr1s', 'nr2s', 'nr3s', 'nspin', 'nqx1', 'nqx2', 'nqx3', 'lda_plus_u_kind', 'edir', 'report', 'esm_nfit', 'space_group', 'origin_choice', 'electron_maxstep', 'mixing_ndim', 'mixing_fixed_ns', 'ortho_para', 'diago_cg_maxiter', 'diago_david_ndim', 'nraise', 'bfgs_ndim', 'if_pos', 'nks', 'nk1', 'nk2', 'nk3', 'sk1', 'sk2', 'sk3', 'nconstr') bool_keys = ('wf_collect', 'tstress', 'tprnfor', 'lkpoint_dir', 'tefield', 'dipfield', 'lelfield', 'lorbm', 'lberry', 'lfcpopt', 'monopole', 'nosym', 'nosym_evc', 'noinv', 'no_t_rev', 'force_symmorphic', 'use_all_frac', 'one_atom_occupations', 'starting_spin_angle', 'noncolin', 'x_gamma_extrapolation', 'lda_plus_u', 'lspinorb', 'london', 'ts_vdw_isolated', 'xdm', 'uniqueb', 'rhombohedral', 'realxz', 'block', 'scf_must_converge', 'adaptive_thr', 'diago_full_acc', 'tqr', 'remove_rigid_rot', 'refold_pos') def smart_int_or_float(numstr): if numstr.find(".") != -1 or numstr.lower().find("e") != -1: return float(numstr) else: return int(numstr) try: if key in bool_keys: if val.lower() == ".true.": return True elif val.lower() == ".false.": return False else: raise ValueError(key + " should be a boolean type!") if key in float_keys: return float(re.search(r"^-?\d*\.?\d*d?-?\d*", val.lower()).group(0).replace("d", "e")) if key in int_keys: return int(re.match(r"^-?[0-9]+", val).group(0)) except ValueError: pass try: val = val.replace("d", "e") return smart_int_or_float(val) except ValueError: pass if "true" in val.lower(): return True if "false" in val.lower(): return False m = re.match(r"^[\"|'](.+)[\"|']$", val) if m: return m.group(1)
[docs]class PWInputError(BaseException): """ Error for PWInput """ pass
[docs]class PWOutput: """ Parser for PWSCF output file. """ patterns = { "energies": r'total energy\s+=\s+([\d\.\-]+)\sRy', "ecut": r'kinetic\-energy cutoff\s+=\s+([\d\.\-]+)\s+Ry', "lattice_type": r'bravais\-lattice index\s+=\s+(\d+)', "celldm1": r"celldm\(1\)=\s+([\d\.]+)\s", "celldm2": r"celldm\(2\)=\s+([\d\.]+)\s", "celldm3": r"celldm\(3\)=\s+([\d\.]+)\s", "celldm4": r"celldm\(4\)=\s+([\d\.]+)\s", "celldm5": r"celldm\(5\)=\s+([\d\.]+)\s", "celldm6": r"celldm\(6\)=\s+([\d\.]+)\s", "nkpts": r"number of k points=\s+([\d]+)" } def __init__(self, filename): """ Args: filename (str): Filename """ self.filename = filename self.data = defaultdict(list) self.read_pattern(PWOutput.patterns) for k, v in self.data.items(): if k == "energies": self.data[k] = [float(i[0][0]) for i in v] elif k in ["lattice_type", "nkpts"]: self.data[k] = int(v[0][0][0]) else: self.data[k] = float(v[0][0][0])
[docs] def read_pattern(self, patterns, reverse=False, terminate_on_match=False, postprocess=str): r""" General pattern reading. Uses monty's regrep method. Takes the same arguments. Args: patterns (dict): A dict of patterns, e.g., {"energy": r"energy\\(sigma->0\\)\\s+=\\s+([\\d\\-.]+)"}. reverse (bool): Read files in reverse. Defaults to false. Useful for large files, esp OUTCARs, especially when used with terminate_on_match. terminate_on_match (bool): Whether to terminate when there is at least one match in each key in pattern. postprocess (callable): A post processing function to convert all matches. Defaults to str, i.e., no change. Renders accessible: Any attribute in patterns. For example, {"energy": r"energy\\(sigma->0\\)\\s+=\\s+([\\d\\-.]+)"} will set the value of self.data["energy"] = [[-1234], [-3453], ...], to the results from regex and postprocess. Note that the returned values are lists of lists, because you can grep multiple items on one line. """ matches = regrep(self.filename, patterns, reverse=reverse, terminate_on_match=terminate_on_match, postprocess=postprocess) self.data.update(matches)
[docs] def get_celldm(self, i): """ Args: i (int): index Returns: Cell dimension along index """ return self.data["celldm%d" % i]
@property def final_energy(self): """ Returns: Final energy """ return self.data["energies"][-1] @property def lattice_type(self): """ Returns: Lattice type. """ return self.data["lattice_type"]