#!/usr/bin/env python
#"""Writen by Niraj K. Nepal, Ph.D."""
"""Module to prepare QE input files"""
import copy
import os
import re
import warnings
from collections import OrderedDict
from pathlib import Path
import numpy as np
import scipy.linalg as alg
from ase.io import espresso
from pymatgen.io.cif import CifWriter
from pymatgen.io.vasp.sets import MPRelaxSet
from pymatgen.io import pwscf
from pymatgen.core.structure import Structure
from pymatgen.core.lattice import Lattice
from pymatgen.core.periodic_table import Element
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from htesp.kpath import kpath
from htesp.check_json import config
from htesp.config import api_key, require_api_key
[docs]
def mprester(key):
"""Return an :class:`MPRester` for ``key``, importing ``mp_api`` on first use.
FIX(20): ``from mp_api.client import MPRester`` at module scope made every
importer of this module -- ``cif_to_gsinput``, ``elastic``, ``site_subs``,
``poscar_to_vasp``, ``qe_input`` -- load ``mp_api``, which in turn imports
``deltalake`` and its compiled extension. Where that extension cannot be
loaded the interpreter *aborts* (``Fatal Python error: Aborted``) rather
than raising, so the whole process dies during the import -- including the
test runner, and including commands that only ever write input files.
Nothing but an actual Materials Project query needs ``mp_api``, so it is
imported here, at the one point where a query is about to happen.
"""
try:
from mp_api.client import MPRester # noqa: PLC0415 -- deliberately lazy
except ImportError as exc: # pragma: no cover - env dependent
raise ImportError(
"Materials Project access needs the mp_api package: "
"pip install mp_api"
) from exc
return MPRester(key)
class _LazyConfig(dict):
"""The configuration, resolved on first use rather than at import time.
``input_data = config()`` at module scope froze the configuration of
whatever directory happened to be current when the module was first
imported. Under the per-material process pool -- and under any
``os.chdir`` -- that is the wrong directory. This proxy resolves (and
re-resolves, when the working directory changes) on the first attribute
access instead, while keeping the ``input_data[...]`` spelling every
function in this module already uses.
"""
_resolved_in = None
def _refresh(self):
here = os.getcwd()
if self._resolved_in != here:
dict.clear(self)
dict.update(self, config())
self._resolved_in = here
return self
def __getitem__(self, key):
return dict.__getitem__(self._refresh(), key)
def get(self, key, default=None):
return dict.get(self._refresh(), key, default)
def __contains__(self, key):
return dict.__contains__(self._refresh(), key)
def keys(self):
return dict.keys(self._refresh())
def items(self):
return dict.items(self._refresh())
input_data = _LazyConfig()
# FIX(5): ``getkpt`` symmetrised with symprec=0.1 and pymatgen's default
# ``international_monoclinic`` while ``setting_qeinput`` used symprec=0.01 with
# ``international_monoclinic=False``, so the k-mesh was computed for a
# different cell than the one written to scf-<mpid>.in. Both go through these
# two constants now.
#: Summary fields that are *nested documents* rather than scalars.
#: FIX(31): ``setting()`` asked for ``available_fields[:-29]``, which still
#: includes these two. emmet-core validates every requested sub-document, and
#: the payload the API returns for them no longer carries the fields its model
#: declares -- ``dos.elemental.Ru.total.1.task_id ... Field required``,
#: ``bandstructure.setyawan_curtarolo.equivalent_labels ... Field required``.
#: One unusable sub-document rejects the *whole* SummaryDoc, so
#: ``mainprogram download`` raised ValidationError for every material that has
#: band-structure or DOS data -- which is nearly all of them -- while the few
#: without it succeeded. Nothing in HTESP reads either field: band structures
#: and densities of states are computed here, not fetched.
UNREQUESTABLE_FIELDS = ("bandstructure", "dos")
#: symmetry tolerance used everywhere a standard cell is derived in this module
SYMPREC = 0.01
#: 2010_SC-style monoclinic setting (alpha < 90, beta = gamma = 90); the
#: international setting puts the unique angle elsewhere, giving a cell that
#: does not match the k-mesh.
INTERNATIONAL_MONOCLINIC = False
# FIX(1c): ``starting_magnetization`` used to be written as 0 for every species
# of an undecorated structure under ``nspin = 2`` -- a non-magnetic calculation
# at twice the cost. This is the last-resort value used when neither the
# structure nor ``config['magmom']['magmom']`` says anything about an element.
DEFAULT_MAGMOM = 0.5
# FIX(2): kept at module level so ``getecut_sssp`` can fall back to it when an
# element is missing from ``config['pseudo']['PSEUDO']``.
#: Wavefunction cutoffs in Ry, mirroring ``pseudo.PSEUDO`` of the packaged
#: config.json so the fallback can never contradict the shipped table (a
#: test asserts they are identical). New elements are taken from SSSP 1.3.0
#: PBE efficiency; the entries that predate it are left as they were.
SSSP_EFFICIENCY = {
'H': 60, 'He': 50, 'Li': 40, 'Be': 40, 'B': 35, 'C': 45, 'N': 60, 'O': 60,
'F': 45, 'Ne': 50, 'Na': 40, 'Mg': 30, 'Al': 30, 'Si': 30, 'P': 30, 'S': 35,
'Cl': 40, 'Ar': 60, 'K': 60, 'Ca': 30, 'Sc': 40, 'Ti': 35, 'V': 35, 'Cr': 40,
'Mn': 65, 'Fe': 90, 'Co': 45, 'Ni': 45, 'Cu': 55, 'Zn': 40, 'Ga': 70, 'Ge': 40,
'As': 35, 'Se': 30, 'Br': 30, 'Kr': 45, 'Rb': 30, 'Sr': 30, 'Y': 35, 'Zr': 30,
'Nb': 40, 'Mo': 35, 'Tc': 30, 'Ru': 35, 'Rh': 35, 'Pd': 45, 'Ag': 50, 'Cd': 60,
'In': 50, 'Sn': 60, 'Sb': 40, 'Te': 30, 'I': 35, 'Xe': 60, 'Cs': 30, 'Ba': 30,
'La': 40, 'Ce': 50, 'Pr': 40, 'Nd': 40, 'Pm': 40, 'Sm': 40, 'Eu': 40, 'Gd': 40,
'Tb': 40, 'Dy': 40, 'Ho': 40, 'Er': 40, 'Tm': 40, 'Yb': 40, 'Lu': 45, 'Hf': 50,
'Ta': 45, 'W': 30, 'Re': 30, 'Os': 40, 'Ir': 55, 'Pt': 35, 'Au': 45, 'Hg': 50,
'Tl': 50, 'Pb': 40, 'Bi': 45, 'Po': 30, 'At': 30, 'Rn': 30, 'Fr': 30, 'Ra': 30,
'Ac': 30, 'Th': 30, 'Pa': 30, 'U': 30, 'Np': 30, 'Pu': 30, 'Am': 30, 'Cm': 30,
'Bk': 60, 'Cf': 60, 'Es': 60, 'Fm': 60, 'Md': 60, 'No': 60, 'Lr': 60,
}
#: start of a QE "card" (everything that is not part of a namelist)
CARD_RE = re.compile(
r"^\s*(ATOMIC_SPECIES|ATOMIC_POSITIONS|ATOMIC_VELOCITIES|ATOMIC_FORCES|"
r"K_POINTS|CELL_PARAMETERS|OCCUPATIONS|CONSTRAINTS|SOLVENTS|HUBBARD)\b")
#: &SYSTEM entries that are rebuilt when a magnetic input is written
SYSTEM_KEY_RE = re.compile(
r"^\s*(nspin|ntyp|nat|starting_magnetization\s*\(\s*\d+\s*\))\s*=")
[docs]
def read_lines(path):
"""Return the lines of ``path`` with their trailing newline stripped."""
return Path(path).read_text().splitlines()
[docs]
def write_lines(path, lines):
"""Write ``lines`` (given without trailing newlines) to ``path``."""
text = "\n".join(lines)
Path(path).write_text(text + "\n" if text else "")
[docs]
def section(lines, start, end):
"""``sed -n '/start/,/end/p' | sed '$d'`` without a shell.
Returns the lines from the first one containing ``start`` up to, but not
including, the first following line containing ``end``. An empty list is
returned when ``start`` never appears.
"""
try:
first = next(i for i, ln in enumerate(lines) if start in ln)
except StopIteration:
return []
for j in range(first + 1, len(lines)):
if end in lines[j]:
return lines[first:j]
return lines[first:]
[docs]
def bare_element(label):
"""``'Fe,spin=5'`` / ``'Fe1'`` -> ``'Fe'``; anything else is returned as is."""
label = str(label).split(",")[0].strip()
match = re.match(r"^([A-Z][a-z]?)", label)
return match.group(1) if match else label
[docs]
def remove_files(*paths):
"""Delete each path that exists, quietly skipping the ones that do not."""
for path in paths:
try:
Path(path).unlink()
except FileNotFoundError:
continue
except OSError as exc:
print("could not remove {}: {}".format(path, exc))
[docs]
def pos_to_kpt(structure_filename,kpoint_density):
"""
Obtain k-point mesh from a structure file.
Parameters:
----------------------
structure_filename : str
Structure file (e.g., QE scf.in or VASP POSCAR).
kpoint_density : float
K-point density.
Returns:
----------------------
kmesh : list
K-point mesh according to the k-point density.
"""
kptsp = kpoint_density
with open(structure_filename,"r") as read_struc:
lines = read_struc.readlines()
# Get cell vectors
line = lines[1].split()
latscl = float(line[0])
ain = np.zeros((3, 3))
amat = np.zeros((3, 3))
bmat = np.zeros((3, 3))
anorm = np.zeros(3)
for i in range(3):
line = lines[2 + i].split()
for j in range(3):
ain[i][j] = float(line[j]) * latscl
amat[j][i] = ain[i][j]
anorm[i] = np.sqrt(ain[i][0] ** 2 + ain[i][1] ** 2 + ain[i][2] ** 2)
bmat = alg.inv(amat)
bnorm = np.zeros(3)
bnorm = alg.norm(bmat,axis=1)
kratio = [bnorm[i] / bnorm[0] for i in range(3)]
klat = bnorm[0] / kptsp
kmesh = [int(kratio[i] * klat + 0.5) if int(kratio[i] * klat + 0.5) != 0 else 1 for i in range(3)]
kratio = [bnorm[i] / bnorm[0] for i in range(3)]
klat = bnorm[0] / kptsp
kmesh = [int(kratio[i] * klat + 0.5) if int(kratio[i] * klat + 0.5) != 0 else 1 for i in range(3)]
return kmesh
[docs]
class MpConnect:
"""
Class for connecting to Materials Project (MP) via MP API, extracting properties,
and preparing input files for Quantum Espresso (QE) calculations.
Parameters:
--------------
key : str, optional
MP API KEY. Use your own key. If not provided, it should be available in config.json
or ../../config.json file.
Attributes:
--------------
prop : list
List of available properties.
mpid : str
Materials ID.
comp : str
Compound name.
data : dict
Dictionary containing materials data.
prefix : str
Prefix for the compound.
ecutwfc : float
Kinetic energy cutoff for wavefunction.
ecutrho : float
Kinetic energy cutoff for charge density.
kpt : list
K-point grid.
structure : Structure
Structure object.
pseudo_dir : str
Path to pseudopotential files.
outdir : str
Directory for output files.
evenkpt : tuple
K-grid with an even number of points in all directions.
kpshift : list
K-point grid shifts.
kptype : str
K-point grid type.
calc : str
Type of calculation.
smear : float
Degauss value for smearing.
smear_type : str
Type of smearing.
etot_conv_thr : float
Total energy convergence threshold.
forc_conv_thr : float
Force convergence threshold.
conv_thr : float
Convergence threshold.
dict_element : dict
Dictionary containing kinetic energy cutoffs for different elements.
comp_list : list
List of elements in the compound.
Example of using MpConnect class:
>>> obj = MpConnect() # Initialize MpConnect object
>>> obj.setting('mp-763') # Set Materials ID to 'mp-763'
>>> obj.getkpt() # Get k-points for the structure
>>> obj.maxecut_sssp() # Calculate maximum recommended energy cutoff using SSSP
>>> obj.setting_qeinput() # Set up Quantum Espresso input files based on the retrieved data
"""
# Intialize the class
def __init__(self):
# The API key comes from $MP_API_KEY, then ~/.config/htesp/credentials,
# then config.json -- never from an unguarded dictionary lookup that
# left ``key`` unbound when no config.json was next to the caller.
# ``api_key`` (not ``require_api_key``) is used here because most users
# of this class -- elastic.py, site_subs.py, poscar_to_vasp.py -- only
# want the input-writing half and never contact Materials Project;
# ``setting`` raises through ``require_api_key`` when a key is needed.
self.key = api_key(input_data)
self.mpr = mprester(self.key) if self.key else None
self.prop = None
self.mpid = None
self.comp = None
self.data = None
self.prefix = None
self.ecutwfc = 0
self.ecutrho = 0
self.kpt = None
self.structure = None
self.pseudo_dir = ""
self.outdir = ""
self.evenkpt = None
self.kpshift = None
self.kptype = ""
self.calc = ""
self.smear = None
self.smear_type = ""
self.etot_conv_thr = None
self.forc_conv_thr = None
self.conv_thr = None
self.dict_element = {}
self.comp_list = None
[docs]
def setting(self,comp):
"""
Initialize the process with a Materials ID.
Parameters:
---------------------
comp : str
Materials ID.
Returns:
---------------------
comp : str
Materials ID.
Notes:
---------------------
This function initializes the process with a Materials ID. It retrieves the data related to the given ID from the Materials Project (MP) database, sets up necessary parameters, and prepares the structure for further calculations.
"""
# FIX(7): kpt, ecutwfc/ecutrho, comp_list and prefix survived from the
# previous material when one process looped over several ids, so a
# material whose own lookup failed silently inherited the last one's
# k-mesh and cutoffs. Reset them before anything else.
self.kpt = None
self.evenkpt = None
self.kptype = ""
self.kpshift = None
self.ecutwfc = 0
self.ecutrho = 0
self.comp_list = None
self.prefix = None
self.structure = None
self.data = None
self.mpid = None
# Set the Material ID
self.comp=comp
if self.mpr is None:
# raises with the instructions for setting $MP_API_KEY
self.key = require_api_key(input_data)
self.mpr = mprester(self.key)
# Retrieve available properties from Materials Project database
# We only store properties except for last 29 elements
# Mostly related to elastic properties, absent for many systems
# FIX(31): drop the nested documents emmet-core cannot validate;
# see UNREQUESTABLE_FIELDS for what that failure looks like.
self.prop = [field
for field in self.mpr.materials.summary.available_fields[:-29]
if field not in UNREQUESTABLE_FIELDS]
# ``config()`` always returns the full schema now, so the settings are
# read unconditionally rather than only when a config.json sits in the
# working directory.
d = input_data["download"]
# Append additional properties specified in the settings to the properties list
for prop in d['element']['prop']:
if prop not in self.prop:
self.prop.append(prop)
#self.data = self.mpr.materials.summary.get_data_by_id(self.comp,fields=self.prop).dict()
# Fetch data for the given Materials ID from the Materials Project database
# FIX(34): ``[0]`` on an empty result is "IndexError: list index out of
# range", which says nothing about which id failed or why. An id that
# Materials Project does not know -- a typo, a deprecated entry, or an
# identifier from another database (examples/QE/tutorial6 ships an
# mpid-list.in full of ``oqmd-*`` ids) -- now names itself.
found = self.mpr.materials.summary.search(material_ids=[self.comp],
fields=self.prop)
if not found:
raise LookupError(
"Materials Project returned no entry for {!r}. Check the id in "
"the tracking file: ids from OQMD ('oqmd-...') or AFLOW "
"('aflow:...') cannot be looked up here, and a deprecated "
"mp-id has to be replaced with its successor.".format(self.comp))
self.data = found[0]
self.structure = self.data.structure
self.data = self.data.dict()
# Extract relevant information from the fetched data
self.mpid = self.data['material_id']
if self.mpid != self.comp:
self.mpid = self.comp
symbol_comp = ""
elm = list(self.data['composition'].keys())
count = list(self.data['composition'].values())
# Construct a prefix based on elemental composition
for i,_ in enumerate(elm):
symbol_comp += elm[i]+str(int(count[i]))
self.prefix = symbol_comp
#print("******************************************\n")
#print("Use property method to extract the following info\n")
#print("property('name') or get_properties(['name1', 'name2', ..])\n")
#print("******************************************\n")
return comp
[docs]
def get_prop_list(self):
"""
Print the list of properties available when downloading the data.
Returns:
---------------------
prop_list : dict_keys
A list of keys representing the available properties.
Notes:
---------------------
This function retrieves the list of properties available for download from the Materials Project (MP) database for the current Materials ID (mpid). It returns a list of keys that represent the available properties that can be downloaded and accessed for analysis or further processing.
"""
#return self.mpr.summary.get_data_by_id(self.mpid).dict().keys()
return self.prop
#return self.mpr.materials.summary.search(material_ids=[self.mpid],fields=self.prop).dict().keys()
[docs]
def download(self,filetype='cif'):
"""
Download the structure file.
Parameters:
---------------
filetype : str, optional
File format for downloading. Default is 'cif'.
Notes:
---------------
This function downloads the structure file for the current Materials ID (mpid) from the Materials Project (MP) database. It saves the file in the specified format, with the default format being CIF ('.cif'). Other supported formats may be specified as needed.
"""
if filetype == 'cif':
unsym_struc = self.structure
CifWriter(unsym_struc, symprec=0.1).write_file('{}.cif'.format(self.mpid))
[docs]
def property(self,name):
"""
Extract a particular property.
Parameters:
------------------
name : str
The name of the property available from the list obtained from the get_prop_list() function.
Returns:
------------------
value
The value of the specified property.
Notes:
------------------
This function extracts the value of a specific property identified by its name from the data retrieved for the current Materials ID (mpid). The property name should be one of the properties listed in the output of the get_prop_list() function.
"""
return self.data[name]
[docs]
def getkpt(self,primitive=True):
"""
Compute k-points based on k-point density.
Parameters:
- primitive (bool): use the primitive standard structure.
Returns:
- kpt (list): the k-point grid.
- kptype (str): the type of k-point grid.
- kptshift (list): the shifts in the k-point grid, returns [0,0,0].
Notes:
---------------------
This function computes the k-point grid based on a specified k-point density,
defaulting to 0.05 if unspecified. It returns the k-point grid, type, and shifts.
"""
if primitive:
# FIX(5): was symprec=0.1 with the international monoclinic
# setting, i.e. a different cell from the one setting_qeinput
# writes; both use the module constants now.
struc = SpacegroupAnalyzer(
self.structure, symprec=SYMPREC
).get_primitive_standard_structure(
international_monoclinic=INTERNATIONAL_MONOCLINIC)
else:
struc = self.structure
relax_set = MPRelaxSet(structure=struc)
relax_set.poscar.write_file('POSCAR')
try:
kptden = input_data['kptden']
except KeyError:
print("Default kptden of 0.05 being utilized\n")
kptden = 0.05
self.kpt = pos_to_kpt("POSCAR",kptden)
self.kptype = "automatic"
self.kpshift = [0, 0, 0]
remove_files("POSCAR")
print("*********************************\n")
print("KPOINT with kpoint density of {} \n".format(kptden))
print("*********************************\n")
return self.kpt,self.kptype,self.kpshift
[docs]
def getevenkpt(self):
"""
Make the k-point mesh even.
Returns:
---------------
evenkpt : tuple
K-grid with an even number of points in all directions.
Notes:
---------------
Ensures the generated k-point grid has even points in each dimension
by incrementing odd components to the next even number.
Returns the resulting even k-point grid as a tuple.
"""
kptsize = len(self.kpt)
kpoint_list = self.kpt
for i in range(kptsize):
if kpoint_list[i]%2 == 0:
kpoint_list[i] = kpoint_list[i]
else:
kpoint_list[i] = kpoint_list[i] + 1
self.evenkpt = tuple(kpoint_list)
return self.evenkpt
[docs]
def getecut_sssp(self,element):
"""
Obtain the kinetic energy cutoff for a specific element.
Parameters:
---------------------
element : str
The type of element for which the kinetic energy cutoff is required.
Returns:
----------------------
float
Kinetic energy cutoff for the particular element.
Raises:
----------------------
KeyError
When the element is in neither ``config['pseudo']['PSEUDO']`` nor
the module's :data:`SSSP_EFFICIENCY` table.
Notes:
----------------------
Gets cutoff values for elements from config.json and falls back to the
SSSP efficiency set for anything the configuration does not cover.
"""
psd_data = input_data.get("pseudo", {}).get("PSEUDO", {})
self.dict_element = dict(psd_data) if psd_data else dict(SSSP_EFFICIENCY)
# spin-decorated species arrive as 'Fe,spin=5'
symbol = bare_element(element)
# FIX(2): the old ``except KeyError`` path re-raised, so one element
# missing from pseudo.PSEUDO aborted the whole scan. Fall back to the
# packaged SSSP table first and only then raise, naming the element.
try:
return self.dict_element[symbol]
except KeyError:
pass
if symbol in SSSP_EFFICIENCY:
warnings.warn(
"element {!r} is missing from pseudo.PSEUDO in config.json; "
"using the SSSP efficiency value of {} Ry".format(
symbol, SSSP_EFFICIENCY[symbol]),
RuntimeWarning, stacklevel=2)
return SSSP_EFFICIENCY[symbol]
raise KeyError(
"no plane-wave cutoff known for element {!r} (from {!r}): it is in "
"neither pseudo.PSEUDO of config.json nor the built-in SSSP "
"efficiency table. Add \"{}\": <ecutwfc in Ry> to the "
"pseudo.PSEUDO dictionary of your config.json.".format(
symbol, element, symbol))
[docs]
def maxecut_sssp(self):
"""
Choose the maximum kinetic energy cutoff among elements in a compound.
Returns:
-----------------
tuple
A tuple containing the maximum kinetic energy cutoff for waveFunction (ecutwfc)
and the maximum kinetic energy cutoff for charge density (ecutrho).
Notes:
-----------------
This function finds the maximum kinetic energy cutoff for compound elements,
then returs the cutoffs for waveFunction and density.
"""
# Try to retrieve elements from the fetched data. The bare ``except:``
# that used to follow swallowed KeyboardInterrupt; ``getecut_sssp``
# now strips any ',spin=' decoration itself, so the third branch is no
# longer needed.
try:
self.comp_list = self.data['elements']
except (KeyError, TypeError):
# If 'elements' is not in the fetched data, take them from the structure
self.comp_list = [str(el) for el in self.structure.elements]
# Get the kinetic energy cutoff for each element in the compound
cutoff_list = [self.getecut_sssp(el) for el in self.comp_list]
# Determine the maximum kinetic energy cutoff for waveFunction
self.ecutwfc= max(cutoff_list)
# Determine the maximum kinetic energy cutoff for charge density
self.ecutrho = 8 * self.ecutwfc
print("******************************************\n")
print("SSSP tested K.E. cutoffs for waveFunction and density\n")
print("******************************************\n")
return self.ecutwfc,self.ecutrho
[docs]
def maxecut_sssp_for_subs(self):
"""
Determine the maximum kinetic energy cutoff among elements in a compound during the substitution process.
Returns:
-----------------
tuple
A tuple containing the maximum kinetic energy cutoff for waveFunction (ecutwfc)
and the maximum kinetic energy cutoff for charge density (ecutrho).
Notes:
-----------------
Similar to maxecut_sssp but for substitution process.
"""
# The bare ``except:`` that used to guard this swallowed
# KeyboardInterrupt and hid a genuinely unknown element;
# ``getecut_sssp`` strips any ',spin=' decoration itself now.
if not self.comp_list:
self.comp_list = [str(el) for el in self.structure.elements]
cutoff_list = [self.getecut_sssp(el) for el in self.comp_list]
self.ecutwfc= max(cutoff_list)
self.ecutrho = 8 * self.ecutwfc
return self.ecutwfc,self.ecutrho
[docs]
def ecut_set(self,ecutwfc=50.0,ecutrho=400.0):
"""
Function to set kinetic energy cutoffs
"""
self.ecutrho=ecutrho
self.ecutwfc=ecutwfc
[docs]
def get_properties(self,property_name):
"""
Function to extract multiple properties from MP.
parameters
-------------------
property_name : list of properties: Default: ['material_id']
Returns
-------------------
property_list : list of properties extracted.
"""
property_list = []
for prop in property_name:
property_list.append(self.data[prop])
property_list.insert(0,self.data['material_id'])
with open('mpid.csv', 'a') as data:
for prop in property_list:
data.write(str(prop) + ",")
data.write("\n")
return property_list
@staticmethod
def _resolve_magnetization(magdict, pseudo_species):
"""Give every *undecorated* species a usable starting magnetization.
Parameters:
- magdict (OrderedDict): label -> sign produced by
:func:`convert_species_list` (0 for a species carrying no spin).
- pseudo_species (list): the original species strings, so a species
that really is spin decorated can be told from one that is not.
Returns:
OrderedDict mapping each label to the value written as
``starting_magnetization``.
"""
# FIX(1c): an undecorated structure got starting_magnetization(i) = 0
# for every species under nspin = 2 -- a non-magnetic calculation at
# twice the cost. Take the per-element value from
# config['magmom']['magmom'] instead, and warn + use DEFAULT_MAGMOM
# when the element is not listed there.
undecorated = {item for item in pseudo_species if ',' not in item}
magmom_cfg = input_data.get('magmom', {}).get('magmom', {}) or {}
resolved = OrderedDict()
for label, value in magdict.items():
if label not in undecorated:
resolved[label] = value
continue
element = bare_element(label)
if element in magmom_cfg:
value = float(magmom_cfg[element])
else:
warnings.warn(
"no spin decoration and no magmom.magmom entry for element "
"{!r}; using starting_magnetization = {} rather than 0, "
"which would make this nspin = 2 run non-magnetic".format(
element, DEFAULT_MAGMOM),
RuntimeWarning, stacklevel=3)
value = DEFAULT_MAGMOM
if abs(value) > 1.0:
warnings.warn(
"magmom.magmom[{!r}] = {} is a moment in Bohr magnetons, but "
"QE's starting_magnetization is a fraction in [-1, 1]; "
"clamping to {}".format(element, value,
1.0 if value > 0 else -1.0),
RuntimeWarning, stacklevel=3)
value = 1.0 if value > 0 else -1.0
resolved[label] = value
return resolved
def _write_magnetic_input(self, scf_name, site_labels, site_coords, magdict):
"""Turn ``scf_name`` into a collinear spin-polarised QE input.
Parameters:
- scf_name (str): the QE input file to rewrite, in place.
- site_labels (list): per-site type label ('Fe1', 'Fe2', 'Pd', ...),
in the site order of the structure that was written.
- site_coords (list): the matching fractional coordinates.
- magdict (dict): label -> starting_magnetization.
The ATOMIC_SPECIES block, ``ntyp``, ``nat``, ``nspin`` and the
``starting_magnetization(i)`` entries are rebuilt from those lists.
"""
lines = read_lines(scf_name)
try:
spc_start = next(i for i, ln in enumerate(lines) if 'ATOMIC_SPECIES' in ln)
pos_start = next(i for i, ln in enumerate(lines) if 'ATOMIC_POSITIONS' in ln)
except StopIteration as exc:
raise RuntimeError(
"{} has no ATOMIC_SPECIES/ATOMIC_POSITIONS block to "
"decorate".format(scf_name)) from exc
natoms = len(site_labels)
if pos_start + natoms >= len(lines):
raise RuntimeError(
"{} has fewer than the {} atomic positions the structure "
"carries".format(scf_name, natoms))
# distinct labels, in order of first appearance -- this is the QE type
# index that starting_magnetization(i) refers to
type_labels = list(OrderedDict.fromkeys(site_labels))
# FIX(1b): one ATOMIC_SPECIES line per *label*, with the pseudopotential
# named after the bare element. The old code renamed line k of the
# block for every spin-decorated species, which only lines up when
# pymatgen happens to emit one line per decorated species, and left
# ntyp describing the undecorated cell.
species_block = ["ATOMIC_SPECIES"]
for label in type_labels:
element = bare_element(label)
species_block.append(" {} {:.4f} {}.upf".format(
label, Element(element).atomic_mass, element))
# FIX(1a): the labels and coordinates come from the standardised cell
position_block = [lines[pos_start]]
for j, label in enumerate(site_labels):
coords = site_coords[j]
position_block.append(" {} {:.10f} {:.10f} {:.10f}".format(
label, coords[0], coords[1], coords[2]))
# end of the existing ATOMIC_SPECIES block
spc_end = spc_start + 1
while (spc_end < len(lines) and lines[spc_end].strip()
and not CARD_RE.match(lines[spc_end])):
spc_end += 1
new_lines = (lines[:spc_start] + species_block
+ lines[spc_end:pos_start] + position_block
+ lines[pos_start + 1 + natoms:])
# drop whatever &SYSTEM already said about spin and type counts
new_lines = [ln for ln in new_lines
if not SYSTEM_KEY_RE.match(ln)]
try:
sys_at = next(i for i, ln in enumerate(new_lines)
if '&SYSTEM' in ln.upper())
except StopIteration as exc:
raise RuntimeError(
"{} has no &SYSTEM namelist".format(scf_name)) from exc
additions = [" nspin = 2,",
" ntyp = {},".format(len(type_labels)),
" nat = {},".format(natoms)]
for index, label in enumerate(type_labels):
additions.append(" starting_magnetization({}) = {},".format(
index + 1, magdict.get(label, DEFAULT_MAGMOM)))
new_lines = new_lines[:sys_at + 1] + additions + new_lines[sys_at + 1:]
write_lines(scf_name, new_lines)
[docs]
def reorder_dictionary(original_dict, new_order):
"""
Reorder the keys of a dictionary based on a new list of keywords and return the list of keys.
Parameters:
- original_dict (dict): The original dictionary whose keys are to be reordered.
- new_order (list): The new list of keywords specifying the desired order of keys.
Returns:
- list: The list of keys of the reordered dictionary based on the new order.
Example:
>>> original_dict = {'a': 1, 'b': 2, 'c': 3, 'd': 4}
>>> new_order = ['c', 'a', 'd', 'b']
>>> reorder_dictionary(original_dict, new_order)
"""
# Create a new OrderedDict with the keys ordered according to the new list
ordered_dict = OrderedDict((key, original_dict[key]) for key in new_order if key in original_dict)
# Return the list of keys in the new order
return ordered_dict
[docs]
def parse_spin(label):
"""The spin out of a species label such as ``'Fe,spin=5'``.
FIX(37): this was ``float(item.split('=')[-1])``, which numpy 2 broke.
``Species.__str__`` renders the spin with ``repr``, and numpy 2 changed the
repr of a scalar from ``5.0`` to ``np.float64(5.0)`` -- so a structure
decorated from ``config['magmom']['magmom']`` (whose values arrive as numpy
floats through pymatgen) produced ``'Fe,spin=np.float64(5.0)'`` and
``mainprogram magenum`` died with
ValueError: could not convert string to float: 'np.float64(5.0)'
Pull the first number out of whatever the label says instead of assuming
the text after ``=`` is already a bare float. Returns 0.0 when the label
carries no number at all.
"""
text = str(label).split("=", 1)[-1].strip()
# 'np.float64(5.0)' -> '5.0'. Taking the first number in the raw text
# would take the 64 out of "float64", so unwrap the call first.
wrapped = re.search(r"\(([^()]*)\)\s*$", text)
if wrapped:
text = wrapped.group(1)
match = re.search(r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?", text)
return float(match.group(0)) if match else 0.0
[docs]
def convert_species_list(original_list):
"""
Rename elements in the original list with elements with different spins.
Parameters:
- original_list (list): List containing strings representing elements.
Returns:
- list: Updated list where elements are replaced with generic names for different spins.
Example:
>>> original_list = ['Fe,spin=5', 'Fe,spin=-5', 'pd', 'pd', 'I,spin=1', 'I,spin=-1']
>>> convert_species_list(original_list)
['Fe1', 'Fe2', 'pd', 'pd', 'I1', 'I2']
"""
elements_spin = extract_elements_with_spin(original_list)
# Updated list
updated_list = []
magdict = {}
# Iterate through the original list
for item in original_list:
if ',' in item:
element = item.split(',')[0]
spin = parse_spin(item)
if len(elements_spin[element]) > 1 and spin > 0:
updated_list.append(element + str(1))
magdict[element + str(1)] = 1
elif len(elements_spin[element]) > 1 and spin < 0:
updated_list.append(element + str(2))
magdict[element + str(2)] = -1
else:
updated_list.append(element)
spin_s = spin
mag_s = 1 if spin_s > 0 else -1 if spin_s < 0 else 0
magdict[element] = mag_s
else:
updated_list.append(item)
magdict[item] = 0
return updated_list,magdict
def ase_cell_to_structure(ase_cell):
# Extract lattice vectors from ASE Atoms object
lattice_vectors = ase_cell.cell.tolist()
# Extract atomic positions and species from ASE Atoms object
atomic_positions = ase_cell.get_scaled_positions()
species_symbols = ase_cell.get_chemical_symbols()
# Create pymatgen Lattice object
lattice = Lattice(lattice_vectors)
# Create list of pymatgen Specie objects
species = [Element(symbol) for symbol in species_symbols]
# Create pymatgen Structure object
structure = Structure(lattice, species, atomic_positions)
return structure
[docs]
def scf_to_dos_scf(input_file='scf.in',out='scf-dos.in'):
"""
Function to change QE input file from using smearing to 'tetrahedra' method, mainly for DOS calculation
parameters
----------------------
input_file : (str) input file for QE scf calculation
out : (str) output file after modification
"""
lines = [ln for ln in read_lines(input_file) if 'calculation' not in ln]
write_lines(input_file, lines)
insert(input_file=input_file,output='temp',keyword='&CONTROL',what="calculation = 'nscf',")
lines = [ln.replace("'smearing'", "'tetrahedra'") for ln in read_lines('temp')]
lines = [ln for ln in lines if 'degaus' not in ln and 'smearing' not in ln]
write_lines(out, lines)
remove_files('temp')
print("Use denser k-mesh for dos calculations")
[docs]
def insert(input_file='scf.in',output='scf-new.in',keyword=None,where="after",what=""):
"""
Function to insert keyword in input file
parameters
-----------------------
input_file : (str) input file. Default: 'scf.in'
output : (str) output file. Default: 'scf-new.in'
keyword : (str) keyword to look after
where : (str) where to insert. Default: 'after', otherwise 'before'
what : (str) what to insert. Default: ''
"""
lines = read_lines(input_file)
out_lines = []
for line in lines:
if keyword is not None and keyword in line:
if where == "after":
out_lines.append(line)
out_lines.append(" " + what)
else:
out_lines.append(" " + what)
out_lines.append(line)
else:
out_lines.append(line)
write_lines(output, out_lines)
# Extract relaxed structure and update QE input.
[docs]
class OUTPUTscf:
"""
class to extract relax structure from QE output file and update input file
parameters
---------------------
filename : (str) QE output file. Default: 'scf.out'
"""
def __init__(self,filename='scf.out'):
self.filename = filename
[docs]
def update_scf(self,input_file='scf.in',output='scf-new.in',what='structure'):
"""
Function to update QE input file and update with new structure
parameters
------------
input_file : (str) QE scf input file to update. Default: 'scf.in'
output : (str) QE scf output file, updated with new structure. Default: 'scf-new.in'
what : (str) what to update. Default: 'structure'. Other, not implemented yet !
"""
lines = read_lines(input_file)
kpoint_block = section(lines, "K_POINTS automatic", "CELL_PARAMETERS angstrom")
header = section(lines, "&CONTROL", "ATOMIC_POSITIONS crystal")
if what == 'structure':
relax = self.extract_relax('relax.in')
write_lines(output, header + kpoint_block + relax)
else:
print('do nothing\n')
remove_files("relax.in")
#def extract_energy(self,input='scf.out'):
# os.system("""Rytoev=`echo "scale=6;13.605698" | bc`""")
# os.system("""en=`grep "! total energy = " {} | tail -n1 | awk '{print $5}'` """.format(input))
# os.system("""e=`echo "scale=6; $en * $Rytoev " | bc`""")
# os.system("""echo "the total energy: $e" """)
# return
#class for calculations. In progress........
#class calculation:
# def __init__(self,batch_header='batch.header'):
# self.batch_header = batch_header
#def create_submission_scripts(self):
#def relax_structure(self,input='scf.in',out='scf.out'):
#def calculate_Tc(self,mu=0.16,degaussq=0.12,ngaussq=0,phfile='elph.out',phsigma=0.005):