Source code for pymatgen.io.pwscf

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

from __future__ import division, unicode_literals

import re
import six

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

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

__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(object): """ 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 == 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 != 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 == None: name = site.specie.symbol+str(c) site_descriptions[name] = site.properties c += 1 def to_str(v): if isinstance(v, six.string_types): 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 == 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) == 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): pass
[docs]class PWOutput(object): 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): 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): """ 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): return self.data["celldm%d" % i]
@property def final_energy(self): return self.data["energies"][-1] @property def lattice_type(self): return self.data["lattice_type"]