"""
Copyright Schrodinger, LLC. All rights reserved.
"""
import argparse
import codecs
from glob import iglob
import json
from itertools import product
import re
import shlex
import warnings
from collections import namedtuple
from html.parser import HTMLParser
from math import cos
from math import log
from math import pi
from math import sin
from math import sqrt
from math import tan
from os import chmod
from os import environ
from os import path
import numpy
from scipy import constants
from pymmlibs import mmfile_get_appdata_release_dir
from schrodinger.application.matsci.nano import space_groups
from schrodinger.application.matsci.nano import xtal
from schrodinger.application.matsci.parserutils import type_positive_int
from schrodinger.application.matsci.msprops import (QE_CART_ATOM_CONSTR,
QE_STARTING_MAG)
from schrodinger.infra import mm
from schrodinger.utils.fileutils import mkdir_p
# File/folder extensions used in IO
# ESM
ESM_EXT = '.esm1'
# PDOS
PROJWFC_UP_EXT = '.projwfc_up'
PROJWFC_DW_EXT = '.projwfc_down'
# NEB
NEB_EXT = '.neb'
# Currently NEB doesn't generate XML, for now copy .out to .nebout
NEB_OUT_EXT = '.nebout'
# SPM file with frequencies and IR intensities (if available)
SPM_EXT = '.vib.spm'
# VIB file with vibrations for Maestro
VIB_EXT = '.vib'
# Force constants file extension (used in matdyn)
FC_EXT = '.fc'
# Frequencies file extension (used in matdyn)
FREQ_EXT = '.freq'
# Phonon band file extension (used in matdyn)
GP_EXT = '%s.gp' % FREQ_EXT
# Phonon DOS file extension (generated by matdyn)
PHDOS_EXT = '.phdos'
# Trajectory directory 'extension'
TRJ_EXT = '_trj'
# Lowdin charges
LOWDIN_EXT = '.lowdin'
# Plotting data (epsilon.x)
DAT_EXT = '.dat'
# epslilon.x prefixes:
EPS_R_PREFIX = 'epsr'
EPS_I_PREFIX = 'epsi'
EPS_PREFIXES = [EPS_R_PREFIX, EPS_I_PREFIX, 'eels', 'ieps', 'jdos']
# turbo_lanczos.x extensions
TURBO_Z_EXTS = tuple('.beta_gamma_z.%d' % idx for idx in range(1, 4))
# turbo_spectrum.x extensions
TURBO_PLOT_EXTS = ('.plot_chi.dat', '.plot_S.dat', '.plot_eps.dat')
# Cube file extension
CUB_EXT = '.cub'
# Hubbard file extension
HP_EXT = '.Hubbard_parameters.dat'
# GIPAW NMR file extension
NMR_EXT = '.nmr.magres'
# Angstrom to Bohr and back
A2B = constants.angstrom / constants.value("Bohr radius")
B2A = 1.0 / A2B
RY2EV = constants.value('Rydberg constant times hc in eV') # From Rydberg to eV
HA2RY = 2.0 # From Hartree to Rydberg
# From THz to cm-1 (approx 33.35641)
THZ2CM1 = constants.value(
'hertz-inverse meter relationship') * constants.tera * constants.centi
# From eV to THz (approx 241.79893)
EV2THZ = constants.value('electron volt-hertz relationship') / constants.tera
# From eV to cm-1 (approx 8065.54429)
EV2CM1 = (constants.value('electron volt-inverse meter relationship') *
constants.centi)
AU2PS = constants.value('atomic unit of time') * 2.0e+12 # From Ry AU to ps
AU2FS = AU2PS * 1.0e+3
# Unit conversion factors: pressure (1 Pa = 1 J/m^3, 1GPa = 10 Kbar )
AU_GPA = (constants.value('Hartree energy') /
(constants.value('Bohr radius')**3 * constants.giga))
HA2KBAR = 10.0 * AU_GPA
PRIMITIVE_PREC = 1e-5
# PDB atom coordinates precision is 1e-3, this should be larger by one order
PDB_PREC = 1e-2
# case insensitive glob filter
UPF_GLOB = '*.[Uu][Pp][Ff]'
UPF_EXT = '.UPF'
# RISM-3D solvent format
MOL_GLOB = '*.MOL'
GREEK_ALPHABET = ('Alpha', 'Beta', 'Gamma', 'Delta', 'Epsilon', 'Zeta', 'Eta',
'Theta', 'Iota', 'Kappa', 'Lambda', 'Mu', 'Nu', 'Xi',
'Omicron', 'Pi', 'Rho', 'Sigma', 'Tau', 'Upsilon')
FLAG_NIMAGE = '-nimage'
FLAG_NPOOLS = '-npools'
FLAG_NTASKS = '-ntg'
FLAG_NBANDS = '-nband'
MagElement = namedtuple('MagElement', ['mag', 'hubb_u', 'hubb_j0', 'zval'])
get_null_mag_element = lambda: MagElement(*([0.] * len(MagElement._fields)))
PseudoData = namedtuple('PseudoData', [
'element', 'ecutwfc', 'ecutrho', 'zval', 'nwfc', 'nbeta', 'pp_type',
'is_fully_rel', 'dft_functional', 'has_gipaw'
],
defaults=['', 0, 0, 0, 0, 0, '', False, '', False])
PPSpecies = namedtuple('PPSpecies', ['ecutwfc', 'ecutrho', 'zval', 'basepath'])
# Regex that matches WFC files
WFC_RE = re.compile(r'(gk*vectors|wfc)[^/]*\.dat$')
# Common command-line flags
FLAG_MPICORES = '-MPICORES'
FLAG_OMP = '-OMP'
SHELL_RUNNER_FILENAME = 'run_qe'
SHELL_RUNNER_DATA = """#!/bin/bash
export TMP_DIR=/scr/$USER
export ESPRESSO_ROOT="$SCHRODINGER"/qe-bin
# Set these variables, if needed
# export PATH=/usr/local/bin:$PATH
# export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
# For OpenMP (serial compilation) set this to empty string
# MPIRUN_EXE=
#
# Use -x flag to export environment variables, if needed
MPIRUN_EXE=orterun
if [ -n "$1" ]; then
exe_name=$1
else
echo "Usage: `basename $0` EXE_NAME TPP INPUT_FILE"
exit 1
fi
if [ -n "$2" ]; then
tpp=$2
# enforce serial execution if variable is -1
if [ "$tpp" -eq "-1" ]; then
tpp="1"
MPIRUN_EXE=
fi
if [ -z "$MPIRUN_EXE" ]; then
export OMP_NUM_THREADS=$tpp
else
export OMP_NUM_THREADS=1
fi
else
echo "Usage: `basename $0` EXE_NAME TPP INPUT_FILE"
exit 1
fi
if [ -n "$3" ]; then
input_file=$3
else
echo "Usage: `basename $0` EXE_NAME TPP INPUT_FILE"
exit 1
fi
EXE_PATH="$ESPRESSO_ROOT/bin/$exe_name"
ARGS=("$@")
P_KEYWORDS=("${ARGS[@]:3}")
echo "===== VARIABLES ====="
echo "EXE_PATH = $EXE_PATH"
echo "tpp = $tpp"
echo "Parallel keywords = ${P_KEYWORDS[@]}"
echo "MPIRUN_EXE = $MPIRUN_EXE"
echo "PATH = $PATH"
echo "LD_LIBRARY_PATH = $LD_LIBRARY_PATH"
if [ -f "$SCHRODINGER_MPI_NODEFILE" ]; then
echo "SCHRODINGER_MPI_NODEFILE = $SCHRODINGER_MPI_NODEFILE:"
cat "$SCHRODINGER_MPI_NODEFILE"
fi
if [ -f "$SCHRODINGER_NODEFILE" ]; then
echo "SCHRODINGER_NODEFILE = $SCHRODINGER_NODEFILE:"
cat "$SCHRODINGER_NODEFILE"
fi
echo "===== END VARIABLES ====="
if [ -z "$MPIRUN_EXE" ]; then
exe_cmd(){ "$EXE_PATH" -in $input_file ; }
else
if [ -f "$SCHRODINGER_MPI_NODEFILE" ]; then
exe_cmd(){ "$MPIRUN_EXE" -hostfile "$SCHRODINGER_MPI_NODEFILE" -np $tpp "$EXE_PATH" -in $input_file "${P_KEYWORDS[@]}" ; }
else
exe_cmd(){ "$MPIRUN_EXE" -np $tpp "$EXE_PATH" -in $input_file "${P_KEYWORDS[@]}" ; }
fi
fi
echo "Invoking: $(declare -f exe_cmd)"
exe_cmd
"""
[docs]def get_shell_runner(check_only=False, write_orig=False, log=None):
"""
Get path to the shell runner script.
:type check_only: bool
:param check_only: If True, only check if runner exists and return
:type write_orig: bool
:param write_orig: If True, create original shell runner even if runner
exists
:type log: function or None
:param log: A function to log info about shell runner existence. Should
take a single str argument.
:rtype: (str, bool)
:return: Path to shell runner and True if shell runner has been just
created, otherwise False
"""
runner_fn = SHELL_RUNNER_FILENAME
status, release_dir = mmfile_get_appdata_release_dir()
mkdir_p(release_dir)
matsci_runner_fn = path.join(release_dir, runner_fn)
matsci_runner_bk_fn = path.join(release_dir, runner_fn + '.orig')
# Paths are sorted by priority
runner_paths = (matsci_runner_fn,
path.join(environ['SCHRODINGER'], 'qe-bin', runner_fn))
write_fns = []
if write_orig:
write_fns += [matsci_runner_bk_fn]
runner_full_path = ''
for runner_fn in runner_paths:
if path.isfile(runner_fn):
runner_full_path = runner_fn
break
if check_only:
if not runner_full_path:
if log:
log('Runner script "run_qe" is not found on this host. Visit '
'https://www.schrodinger.com/kb/1954 for installation '
'instructions.\nTo generate "run_qe", run:\n'
'$SCHRODINGER/run periodic_dft_gui_dir/runner.py -h')
return runner_full_path, False
if not runner_full_path:
runner_full_path = matsci_runner_fn
write_fns += [runner_full_path]
for runner_fn in write_fns:
with open(runner_fn, 'w') as runner_fh:
runner_fh.write(SHELL_RUNNER_DATA)
if not runner_fn.endswith('.bk'):
chmod(runner_fn, 0o755)
runner_written = len(write_fns) == 2
return runner_full_path, runner_written
[docs]class ArgumentParserNoExit(argparse.ArgumentParser):
"""
Subclass that raises instead of exiting on error.
"""
[docs] def error(self, message):
"""
From python docs: If you override this in a subclass, it should not
return -- it should either exit or raise an exception.
"""
raise argparse.ArgumentTypeError(message)
[docs]def add_qe_parallel_parser_arguments(parser):
"""
Add QE parallel arguments to the parser.
:type parser: `argparse.ArgumentParser`
:param parser: The parser to add arguments to
"""
rsopts = parser.add_argument_group('Quantum ESPRESSO parallel options')
rsopts.add_argument(
FLAG_NIMAGE,
type=type_positive_int,
default=1,
help='Number of images (points in the configuration space)')
rsopts.add_argument(FLAG_NPOOLS,
type=type_positive_int,
default=1,
help='Number of K-point groups to run in parallel')
group = rsopts.add_mutually_exclusive_group()
group.add_argument(FLAG_NTASKS,
type=type_positive_int,
help='Number of task groups to run in parallel')
group.add_argument(FLAG_NBANDS,
type=type_positive_int,
help='Number of bands to run in parallel')
[docs]def get_pkeywords(options, ncpus, npools=None):
"""
Validate and get string of parallel keywords for QE binaries (pw.x, etc.).
:type options: argparse Namespace object
:param options: The input options
:type ncpus: int
:param ncpus: Number of total requested CPUs
:param int npools: npools value, None for default passed in options
:rtype: bool, str or list
:return: If converted successfully, True and a list of [keyword, value] are
returned, otherwise, False and error message are returned
"""
msg = '%s value (%d) is not a divisor of number of processors (%d)'
npools = npools if npools is not None else options.npools
if ncpus % options.nimage != 0:
return (False, msg % (FLAG_NIMAGE, options.nimage, ncpus))
# Get quotient
ncpus_per_image = ncpus // options.nimage
if ncpus_per_image % npools != 0:
return (False, msg % (FLAG_NPOOLS, npools, ncpus_per_image))
cmd = [FLAG_NIMAGE, str(options.nimage), FLAG_NPOOLS, str(npools)]
if options.ntg:
ncpus_per_pool = ncpus_per_image // npools
if ncpus_per_pool % options.ntg != 0:
return (False, msg % (FLAG_NTASKS, options.ntg, ncpus_per_pool))
cmd += [FLAG_NTASKS, str(options.ntg)]
elif options.nband:
ncpus_per_pool = ncpus_per_image // npools
if ncpus_per_pool % options.nband != 0:
return (False, msg % (FLAG_NBANDS, options.nband, ncpus_per_pool))
cmd += [FLAG_NBANDS, str(options.nband)]
return True, cmd
[docs]def get_mag_hubbu(atom, decimals=10):
"""
Get starting magnetization and Hubbard U parameters from atom property. If
not present, return the default.
:type atom: `schrodinger.structure._StructureAtom`
:param atom: Atom from which values are taken
:type decimals: int
:param decimals: Number of decimals to round to
:rtype: MagElement namedtuple
:return: Tuple of starting magnetization and Hubbard parameters
"""
mag_zero = get_null_mag_element()
# round is important since structure is a C object and floats might be
# represented weird in C vs Python
mag_element = numpy.round(
json.loads(atom.property.get(QE_STARTING_MAG, json.dumps(mag_zero))),
decimals).tolist()
assert len(mag_element) <= len(MagElement._fields)
if len(mag_element) < len(MagElement._fields):
# Padding with zeros for backward compatibility
mag_element.extend([0] * (len(MagElement._fields) - len(mag_element)))
return MagElement(*mag_element)
[docs]def has_hubbard(struct):
"""
Check if input structure has at least one atom with Hubbard U set.
:param structure.Structure struct: Input structure
:return bool: Whether at least one atom has Hubbard U set
"""
for atom in struct.atom:
mag_hubbu = get_mag_hubbu(atom)
if mag_hubbu.hubb_u or mag_hubbu.hubb_j0:
return True
return False
[docs]def set_mag_hubbu(atom, mag_element):
"""
Set starting magnetization and Hubbard U parameters in atom property.
:type atom: `schrodinger.structure._StructureAtom`
:param atom: Atom for which values are set
:type mag_element: MagElement namedtuple
:param mag_element: Tuple of starting magnetization and Hubbard parameters
"""
# Convert namedtuple to list for json
atom.property[QE_STARTING_MAG] = json.dumps(list(mag_element))
[docs]def sync_pbc(struct,
lattice_params=None,
chorus_params=None,
prioritize_cparams=False):
"""
Sync PBC properties in place (without creating a new structure) for struct.
If all PBC properties are absent (both chorus and PDB) return False. It is
possible to provide new lattice or chorus parameters and to prioritize one
of the sets.
:type: `schrodinger.structure.Structure`
:param: Structure to be modified
:type lattice_params: list or numpy.array or None
:param lattice_params: contains the a, b, c, alpha, beta, and gamma lattice
parameters or None. These will be used instead of ones possibly obtained
from the struct
:type chorus_params: list or numpy.array or None
:param chorus_params: contains the nine chorus properties,
i.e. ax, ay, az, bx, ..., cz or None. These will be used instead of ones
possibly obtained from the struct
:type prioritize_cparams: bool
:param: Prioritize chorus params over lattice params if True. If False,
lattice params have priority
:rtype: bool
:return: True on syncing success, False if both sets were not provided or
empty
"""
try:
if lattice_params is None:
lattice_params = xtal.get_lattice_param_properties(struct)
except KeyError:
lattice_params = []
try:
if chorus_params is None:
chorus_params = xtal.get_chorus_properties(struct)
except KeyError:
chorus_params = []
if not any(map(len, [lattice_params, chorus_params])):
return False
use_cparams = prioritize_cparams and len(chorus_params)
if not use_cparams and len(lattice_params):
xtal.set_pbc_properties(struct,
xtal.get_chorus_from_params(lattice_params))
else:
xtal.set_pbc_properties(struct, chorus_params)
return True
[docs]def get_mpircores_from_environ():
"""
Get -MPICORES flag value from SCHRODINGER_COMMANDLINE environment variable.
If not found, return 1.
:rtype: int
:return: MPICORES value or 1 (if not defined)
"""
cmds = shlex.split(environ['SCHRODINGER_COMMANDLINE'])
for idx, cmd in enumerate(cmds):
if cmd == '-MPICORES' and len(cmds) > idx + 1:
return int(cmds[idx + 1])
return 1
[docs]def get_element(element):
"""
Get element name with the removed digits from atom type. Those might be
present to due starting magnetization.
:type element: str
:param element: Element
:rtype: str
:return: Element without digits
"""
return ''.join(x for x in element if not x.isdigit())
[docs]def copy_atom_constr(src_atom, dest_atom):
"""
Copy Cartesian atomic constraints (if any) from source atom to destination
atom.
:param structure._StructureAtom src_atom: Source atom
:param structure._StructureAtom dest_atom: Destination atom
"""
val = src_atom.property.get(QE_CART_ATOM_CONSTR)
if val is not None:
dest_atom.property[QE_CART_ATOM_CONSTR] = val
[docs]def set_atom_cart_constr(atom, constr):
"""
"""
if isinstance(constr, str):
atom.property[QE_CART_ATOM_CONSTR] = constr
return
constr = list(map(int, constr))
assert len(constr) == 3
if not all(constr):
atom.property[QE_CART_ATOM_CONSTR] = json.dumps(constr)
[docs]def set_job_builder_opts(job_builder, cfg_fn=None, add_pps=False):
"""
Set options to job builder.
:param launchapi.JobSpecificationArgsBuilder job_builder: Job builder to
update
:param str cfg_fn: Config file to set as input
:param bool add_pps: Whether to add pseudopotential files
"""
if cfg_fn:
job_builder.setInputFile(cfg_fn)
if add_pps:
for upf_fn in iglob(UPF_GLOB):
job_builder.setInputFile(upf_fn)
[docs]class UPFParser(object):
"""
Class that handles UPF parsing.
"""
PP_TYPE_NC, PP_TYPE_PAW = 'NC', 'PAW'
PP_TYPES = ('1/r', 'US', PP_TYPE_NC, PP_TYPE_PAW)
FULLY_REL = 'full'
UPF2_TRUE = ['T']
[docs] def __init__(self, path, file_fh=None, binary=False):
"""
Initialize UPFParser.
:type path: str
:param path: Path to the UPF file
"""
assert any([path, file_fh])
self.pseudo = PseudoData()
try:
if path:
upf_fh = open(path)
else:
upf_fh = file_fh
if binary:
reader = codecs.getreader('utf-8')
upf_fh = reader(upf_fh)
first_line = upf_fh.readline().strip()
if '<?xml' in first_line:
# Skip XML header, might be present or not
first_line = upf_fh.readline().strip()
if ('<UPF version="2.0.0">' in first_line or
'<UPF version="2.0.1">' in first_line):
self.pseudo = self._parseUPF2(upf_fh)
elif '<PP_INFO>' in first_line:
self.pseudo = self._parseUPF(upf_fh)
except (IOError, TypeError, IndexError, ValueError):
pass
finally:
if path:
upf_fh.close()
[docs] def getPseudo(self):
""" Return a copy of the pseudo data."""
return PseudoData._make(self.pseudo)
@staticmethod
def _parseUPF2Prop(line):
"""
Get value of the property from line of a UPF2 file.
:type line: str
:param: Line from UPF2 file with a property
:rtype: str
:return: Parsed property
"""
return line.split('=')[-1].strip().strip('"\' ')
@staticmethod
def _parseUPF2(upf_fh):
"""
Parse UPF version 2 file and set related attributes.
:type upf_fh: file handler
:param upf_fh: File handler to UPF
"""
pp_header_parser = PPHeaderHTMLParser()
in_pp_header = False
for line in upf_fh:
if '<PP_MESH' in line:
pp_header_parser.close()
break
if '<PP_HEADER' in line:
pp_header_parser.feed(line)
in_pp_header = True
continue
if in_pp_header:
pp_header_parser.feed(line)
data = dict(PseudoData._field_defaults)
if 'wfc_cutoff' in pp_header_parser.data:
data['ecutwfc'] = float(pp_header_parser.data['wfc_cutoff'])
if 'rho_cutoff' in pp_header_parser.data:
data['ecutrho'] = float(pp_header_parser.data['rho_cutoff'])
if 'z_valence' in pp_header_parser.data:
data['zval'] = float(pp_header_parser.data['z_valence'])
if 'element' in pp_header_parser.data:
element = pp_header_parser.data['element'].title()
data['element'] = UPFParser._checkElement(element)
if 'pseudo_type' in pp_header_parser.data:
pp_type = pp_header_parser.data['pseudo_type']
data['pp_type'] = UPFParser._checkType(pp_type)
if 'relativistic' in pp_header_parser.data:
is_fully_rel = pp_header_parser.data['relativistic']
data['is_fully_rel'] = UPFParser._checkRelativistic(is_fully_rel)
if 'functional' in pp_header_parser.data:
dft_functional = pp_header_parser.data['functional']
data['dft_functional'] = UPFParser._getFunctional(dft_functional)
if 'has_gipaw' in pp_header_parser.data:
data['has_gipaw'] = (pp_header_parser.data['has_gipaw']
in UPFParser.UPF2_TRUE)
return PseudoData(**data)
@staticmethod
def _parseUPFProp(line, get_first=True):
"""
Get value of the property from line of a UPF file.
:type line: str
:param line: Line from UPF file with a property
:type get_first: bool
:param get_first: If True, strip and split line, otherwise only strip
:rtype: str
:return: Parsed property
"""
ret = line.strip()
if get_first:
ret = ret.split()[0].strip()
return ret
@staticmethod
def _parseUPFAddInfo(nwfc, nbeta, upf_fh):
"""
Parse additional info section of the UPF v1 file. The file handler is
expected to be at the '<PP_ADDINFO>' line. Sets self.is_fully_rel to
True, if pseudopotential is fully relativistic.
:type upf_fh: file handler
:param upf_fh: File handler to UPF
"""
# It is not fully relativistic if j != l+1/2 and j != l-1/2
for idx in range(nwfc):
tmp = upf_fh.readline().strip().split()
lchi, jchi = list(map(float, tmp[2:4]))
diff = jchi - lchi
if abs(diff - 0.5) > 1e-7 and abs(diff + 0.5) > 1e-7:
return False
for idx in range(nbeta):
tmp = upf_fh.readline().strip().split()
lll, jjj = list(map(float, tmp[0:2]))
diff = lll - jjj
if abs(diff - 0.5) > 1e-7 and abs(diff + 0.5) > 1e-7:
return False
return True
@staticmethod
def _parseUPF(upf_fh):
"""
Parse UPF version 1 file and set related attributes.
:type upf_fh: file handler
:param upf_fh: File handler to UPF
"""
data = dict(PseudoData._field_defaults)
while True:
line = upf_fh.readline()
if not line:
break
if '<PP_ADDINFO' in line:
data['is_fully_rel'] = UPFParser._parseUPFAddInfo(
data['nwfc'], data['nbeta'], upf_fh)
if '<PP_GIPAW_FORMAT_VERSION' in line:
data['has_gipaw'] = True
if '<PP_HEADER' in line:
# Start reading UPF v1 parameters
# Skip version number
upf_fh.readline()
element = UPFParser._parseUPFProp(upf_fh.readline()).title()
data['element'] = UPFParser._checkElement(element)
# Pseudopotential type
pp_type = UPFParser._parseUPFProp(upf_fh.readline())
data['pp_type'] = UPFParser._checkType(pp_type)
# Skip NLCC
upf_fh.readline()
# DFT functional
dft_functional = UPFParser._parseUPFProp(upf_fh.readline(),
get_first=False)
data['dft_functional'] = UPFParser._getFunctional(
dft_functional)
# Z val
data['zval'] = float(UPFParser._parseUPFProp(upf_fh.readline()))
# Skip total energy
upf_fh.readline()
# Ecutwfc and Ecutrho
tmp = upf_fh.readline().strip().split()
data['ecutwfc'] = float(tmp[0].strip())
data['ecutrho'] = float(tmp[1].strip())
# Skip angular momentum
upf_fh.readline()
# Skip mumber of points in mesh
upf_fh.readline()
# Number of Wavefunctions, Number of Projectors
tmp = upf_fh.readline().strip().split()
data['nwfc'], data['nbeta'] = list(map(int, tmp[:2]))
return PseudoData(**data)
@staticmethod
def _checkElement(element):
"""
Check that element is known by mm infrastructure.
:type element: str
:param element: Element name
:rtype: str
:return: If element is know return element name, otherwise empty string
"""
atomic_number = mm.mmelement_get_atomic_number_by_symbol(element)
if atomic_number > 0:
return element
else:
return ''
@staticmethod
def _checkType(pp_type):
"""
Check that pseudopotential type is known by our infrastructure.
:type pp_type: str
:param pp_type: Pseudopotential type
:rtype: str
:return: If pseudopotential type is know return pseudopotential type,
otherwise empty string
"""
for pp_allowed_type in UPFParser.PP_TYPES:
# This deals with 'USPP' instead of just 'US'
if pp_type.upper().startswith(pp_allowed_type):
return pp_allowed_type
# 'SL' is the same as 'NC'
elif pp_type.upper().startswith('SL'):
return 'NC'
return ''
@staticmethod
def _checkRelativistic(rel_str):
"""
Check if pseudo is fully relativistic.
:type rel_str: str
:param rel_str: Relativistic type
:rtype: bool
:return: If pseudo is fully relativistic return True, otherwise False
"""
return rel_str == UPFParser.FULLY_REL
@staticmethod
def _getFunctional(func_str):
"""
Get DFT functional based on the string from PP file
:type func_str: str
:param func_str: String describing DFT functional
:rtype: str
:return: DFT functional
"""
if any(func in func_str
for func in ['PBESOL', 'SLA-PW-PSX-PSC', 'SLA PW PSX PSC']):
return 'PBESOL'
elif any(func in func_str
for func in ['PBE', 'SLA-PW-PBX-PBC', 'SLA PW PBX PBC']):
return 'PBE'
elif any(
func in func_str
for func in ['LDA', 'SLA-PZ-NOGX-NOGC', 'SLA PZ NOGX NOGC']):
return 'PZ'
elif any(func in func_str
for func in ['PW91', 'SLA-PW-GGX-GGC', 'SLA PW GGX GGC']):
return 'PW91'
elif 'BP' in func_str:
return 'BP'
elif 'BLYP' in func_str:
return 'BLYP'
return ''
[docs]class HighSymmetryKPath:
"""
Based on symmetry/bandstructure.py :: HighSymmKpath (MIT license)
This class looks for path along high symmetry lines in
the Brillouin Zone.
It is based on Setyawan, W., & Curtarolo, S. (2010).
High-throughput electronic band structure calculations:
Challenges and tools. Computational Materials Science,
49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010
"""
TWOD = '_2D'
[docs] def __init__(self, struct, is_2d=False):
"""
Initialize kpath class. self.kpath will contain 'edge' points.
:param structure.Structure struct: Structure that has space group or
unit cell data
:param bool is_2d: Whether to build K-point path for 2D/ESM calcs
"""
self.kpath = []
# Let it raise KeyError, if necessary
spg_symbol = struct.property[xtal.SPACE_GROUP_KEY]
spg_obj = space_groups.get_spacegroups().getSpgObjByName(spg_symbol)
if not spg_obj:
raise ValueError("Couldn't find lattice type for the %s symbol" %
spg_symbol)
# Always use short name
spg_symbol = spg_obj.space_group_short_name
latt_type = spg_obj.crystal_system.name
cell_dimensions = xtal.get_lattice_param_properties(struct)
a_latt, b_latt, c_latt, alpha, beta, gamma = cell_dimensions
alpha_r, beta_r, gamma_r = numpy.radians(cell_dimensions[3:])
error_msg = ('Unexpected space group symbol: %s, for %s lattice type.' %
(spg_symbol, latt_type))
if latt_type == space_groups.CrystalSystems.CUBIC_NAME:
if 'P' in spg_symbol:
kpath = self.cubic()
elif 'F' in spg_symbol:
kpath = self.fcc()
elif 'I' in spg_symbol:
kpath = self.bcc()
else:
raise ValueError(error_msg)
elif latt_type == space_groups.CrystalSystems.TETRAGONAL_NAME:
if 'P' in spg_symbol:
kpath = self.tet()
elif 'I' in spg_symbol:
if c_latt < a_latt:
kpath = self.bctet1(c_latt, a_latt)
else:
kpath = self.bctet2(c_latt, a_latt)
else:
raise ValueError(error_msg)
elif latt_type == space_groups.CrystalSystems.ORTHORHOMBIC_NAME:
if 'F' in spg_symbol:
if 1 / a_latt**2 > 1 / b_latt**2 + 1 / c_latt**2:
kpath = self.orcf1(a_latt, b_latt, c_latt)
elif 1 / a_latt**2 < 1 / b_latt**2 + 1 / c_latt**2:
kpath = self.orcf2(a_latt, b_latt, c_latt)
else:
kpath = self.orcf3(a_latt, b_latt, c_latt)
elif 'I' in spg_symbol:
kpath = self.orci(a_latt, b_latt, c_latt)
elif 'C' in spg_symbol:
kpath = self.orcc(a_latt, b_latt, c_latt)
else:
kpath = self.orc()
elif latt_type == space_groups.CrystalSystems.HEXAGONAL_NAME:
kpath = self.hex(is_2d)
elif latt_type == space_groups.CrystalSystems.TRIGONAL_NAME:
if spg_symbol.startswith('P'):
# Trigonal P space groups are hexagonal
kpath = self.hex(is_2d)
elif alpha < 90.0:
kpath = self.rhl1(alpha_r)
else:
kpath = self.rhl2(alpha_r)
elif latt_type == space_groups.CrystalSystems.MONOCLINIC_NAME:
if 'C' in spg_symbol:
kpath = self.mclc1(a_latt, b_latt, c_latt, alpha_r)
else:
kpath = self.mcl(b_latt, c_latt, alpha_r)
# Triclinic cell remained
else:
kpath = self.tria(is_2d)
if is_2d and not self.name.endswith(self.TWOD):
# Enforce 2D path if requested
kpath = self.tria(is_2d)
for name in kpath['path']:
coords = kpath['kpoints'][name]
self.kpath.append(coords.tolist() + [name])
[docs] def cubic(self):
"""
Generate 'edge' k-points for the cubic lattice.
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "CUB"
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'X': numpy.array([0.0, 0.5, 0.0]),
'R': numpy.array([0.5, 0.5, 0.5]),
'M': numpy.array([0.5, 0.5, 0.0])
}
path = [r"\Gamma", "X", "M", r"\Gamma", "R", "X", "M", "R"]
return {'kpoints': kpoints, 'path': path}
[docs] def fcc(self):
"""
Generate 'edge' k-points for the FCC cubic lattice.
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "FCC"
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'K': numpy.array([0.375, 0.375, 0.75]),
'L': numpy.array([0.5, 0.5, 0.5]),
'U': numpy.array([0.625, 0.25, 0.625]),
'W': numpy.array([0.5, 0.25, 0.75]),
'X': numpy.array([0.5, 0.0, 0.5])
}
path = [
r"\Gamma", "X", "W", "K", r"\Gamma", "L", "U", "W", "L", "K", "U",
"X"
]
return {'kpoints': kpoints, 'path': path}
[docs] def bcc(self):
"""
Generate 'edge' k-points for the BCC cubic lattice.
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "BCC"
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'H': numpy.array([0.5, -0.5, 0.5]),
'P': numpy.array([0.25, 0.25, 0.25]),
'N': numpy.array([0.0, 0.0, 0.5])
}
path = [r"\Gamma", "H", "N", r"\Gamma", "P", "H", "P", "N"]
return {'kpoints': kpoints, 'path': path}
[docs] def tet(self):
"""
Generate 'edge' k-points for the tetrahedral lattice.
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "TET"
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'A': numpy.array([0.5, 0.5, 0.5]),
'M': numpy.array([0.5, 0.5, 0.0]),
'R': numpy.array([0.0, 0.5, 0.5]),
'X': numpy.array([0.0, 0.5, 0.0]),
'Z': numpy.array([0.0, 0.0, 0.5])
}
path = [
r"\Gamma", "X", "M", r"\Gamma", "Z", "R", "A", "Z", "X", "R", "M",
"A"
]
return {'kpoints': kpoints, 'path': path}
[docs] def bctet1(self, c, a):
"""
Generate 'edge' k-points for the tetrahedral lattice with I setting.
:type c: float
:param c: C length
:type a: float
:param a: A length
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "BCT1"
eta = (1 + c**2 / a**2) / 4
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'M': numpy.array([-0.5, 0.5, 0.5]),
'N': numpy.array([0.0, 0.5, 0.0]),
'P': numpy.array([0.25, 0.25, 0.25]),
'X': numpy.array([0.0, 0.0, 0.5]),
'Z': numpy.array([eta, eta, -eta]),
'Z_1': numpy.array([-eta, 1 - eta, eta])
}
path = [
r"\Gamma", "X", "M", r"\Gamma", "Z", "P", "N", "Z_1", "M", "X", "P"
]
return {'kpoints': kpoints, 'path': path}
[docs] def bctet2(self, c, a):
"""
Generate 'edge' k-points for the tetrahedral lattice with I setting.
:type c: float
:param c: C length
:type a: float
:param a: A length
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "BCT2"
eta = (1 + a**2 / c**2) / 4.0
zeta = a**2 / (2.0 * c**2)
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'N': numpy.array([0.0, 0.5, 0.0]),
'P': numpy.array([0.25, 0.25, 0.25]),
r'\Sigma': numpy.array([-eta, eta, eta]),
r'\Sigma_1': numpy.array([eta, 1 - eta, -eta]),
'X': numpy.array([0.0, 0.0, 0.5]),
'Y': numpy.array([-zeta, zeta, 0.5]),
'Y_1': numpy.array([0.5, 0.5, -zeta]),
'Z': numpy.array([0.5, 0.5, -0.5])
}
path = [
r"\Gamma", "X", "Y", r"\Sigma", r"\Gamma", "Z", r"\Sigma_1", "N",
"P", "Y_1", "Z", "X", "P"
]
return {'kpoints': kpoints, 'path': path}
[docs] def orc(self):
"""
Generate 'edge' k-points for the orthorhombic lattice.
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "ORC"
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'R': numpy.array([0.5, 0.5, 0.5]),
'S': numpy.array([0.5, 0.5, 0.0]),
'T': numpy.array([0.0, 0.5, 0.5]),
'U': numpy.array([0.5, 0.0, 0.5]),
'X': numpy.array([0.5, 0.0, 0.0]),
'Y': numpy.array([0.0, 0.5, 0.0]),
'Z': numpy.array([0.0, 0.0, 0.5])
}
path = [
r"\Gamma", "X", "S", "Y", r"\Gamma", "Z", "U", "R", "T", "Z", "Y",
"T", "U", "X", "S", "R"
]
return {'kpoints': kpoints, 'path': path}
[docs] def orcf1(self, a, b, c):
"""
Generate 'edge' k-points for the orthorhombic lattice with F setting.
:type a: float
:param a: A length
:type b: float
:param b: B length
:type c: float
:param c: C length
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "ORCF1"
zeta = (1 + a**2 / b**2 - a**2 / c**2) / 4.0
eta = (1 + a**2 / b**2 + a**2 / c**2) / 4.0
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'A': numpy.array([0.5, 0.5 + zeta, zeta]),
'A_1': numpy.array([0.5, 0.5 - zeta, 1 - zeta]),
'L': numpy.array([0.5, 0.5, 0.5]),
'T': numpy.array([1.0, 0.5, 0.5]),
'X': numpy.array([0.0, eta, eta]),
'X_1': numpy.array([1, 1 - eta, 1 - eta]),
'Y': numpy.array([0.5, 0.0, 0.5]),
'Z': numpy.array([0.5, 0.5, 0.0])
}
path = [
r"\Gamma", "Y", "T", "Z", r"\Gamma", "X", "A_1", "Y", "T", "X_1",
"X", "A", "Z", "L", r"\Gamma"
]
return {'kpoints': kpoints, 'path': path}
[docs] def orcf2(self, a, b, c):
"""
Generate 'edge' k-points for the orthorhombic lattice with F setting.
:type a: float
:param a: A length
:type b: float
:param b: B length
:type c: float
:param c: C length
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "ORCF2"
phi = (1 + c**2 / b**2 - c**2 / a**2) / 4.0
eta = (1 + a**2 / b**2 - a**2 / c**2) / 4.0
delta = (1 + b**2 / a**2 - b**2 / c**2) / 4.0
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'C': numpy.array([0.5, 0.5 - eta, 1 - eta]),
'C_1': numpy.array([0.5, 0.5 + eta, eta]),
'D': numpy.array([0.5 - delta, 0.5, 1 - delta]),
'D_1': numpy.array([0.5 + delta, 0.5, delta]),
'L': numpy.array([0.5, 0.5, 0.5]),
'H': numpy.array([1 - phi, 0.5 - phi, 0.5]),
'H_1': numpy.array([phi, 0.5 + phi, 0.5]),
'X': numpy.array([0.0, 0.5, 0.5]),
'Y': numpy.array([0.5, 0.0, 0.5]),
'Z': numpy.array([0.5, 0.5, 0.0])
}
path = [
r"\Gamma", "Y", "C", "D", "X", r"\Gamma", "Z", "D_1", "H", "C",
"C_1", "Z", "X", "H_1", "H", "Y", "L", r"\Gamma"
]
return {'kpoints': kpoints, 'path': path}
[docs] def orcf3(self, a, b, c):
"""
Generate 'edge' k-points for the orthorhombic lattice with F setting.
:type a: float
:param a: A length
:type b: float
:param b: B length
:type c: float
:param c: C length
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "ORCF3"
zeta = (1 + a**2 / b**2 - a**2 / c**2) / 4.0
eta = (1 + a**2 / b**2 + a**2 / c**2) / 4.0
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'A': numpy.array([0.5, 0.5 + zeta, zeta]),
'A_1': numpy.array([0.5, 0.5 - zeta, 1 - zeta]),
'L': numpy.array([0.5, 0.5, 0.5]),
'T': numpy.array([1.0, 0.5, 0.5]),
'X': numpy.array([0.0, eta, eta]),
'X_1': numpy.array([1, 1 - eta, 1 - eta]),
'Y': numpy.array([0.5, 0.0, 0.5]),
'Z': numpy.array([0.5, 0.5, 0.0])
}
path = [
r"\Gamma", "Y", "T", "Z", r"\Gamma", "X", "A_1", "Y", "X", "A", "Z",
"L", r"\Gamma"
]
return {'kpoints': kpoints, 'path': path}
[docs] def orci(self, a, b, c):
"""
Generate 'edge' k-points for the orthorhombic lattice with I setting.
:type a: float
:param a: A length
:type b: float
:param b: B length
:type c: float
:param c: C length
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "ORCI"
zeta = (1 + a**2 / c**2) / 4.0
eta = (1 + b**2 / c**2) / 4.0
delta = (b**2 - a**2) / (4 * c**2)
mu = (a**2 + b**2) / (4 * c**2)
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'L': numpy.array([-mu, mu, 0.5 - delta]),
'L_1': numpy.array([mu, -mu, 0.5 + delta]),
'L_2': numpy.array([0.5 - delta, 0.5 + delta, -mu]),
'R': numpy.array([0.0, 0.5, 0.0]),
'S': numpy.array([0.5, 0.0, 0.0]),
'T': numpy.array([0.0, 0.0, 0.5]),
'W': numpy.array([0.25, 0.25, 0.25]),
'X': numpy.array([-zeta, zeta, zeta]),
'X_1': numpy.array([zeta, 1 - zeta, -zeta]),
'Y': numpy.array([eta, -eta, eta]),
'Y_1': numpy.array([1 - eta, eta, -eta]),
'Z': numpy.array([0.5, 0.5, -0.5])
}
path = [
r"\Gamma", "X", "L", "T", "W", "R", "X_1", "Z", r"\Gamma", "Y", "S",
"W", "L_1", "Y", "Y_1", "Z"
]
return {'kpoints': kpoints, 'path': path}
[docs] def orcc(self, a, b, c):
"""
Generate 'edge' k-points for the orthorhombic lattice with C setting.
:type a: float
:param a: A length
:type b: float
:param b: B length
:type c: float
:param c: C length
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "ORCC"
zeta = (1 + a**2 / b**2) / 4.0
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'A': numpy.array([zeta, zeta, 0.5]),
'A_1': numpy.array([-zeta, 1 - zeta, 0.5]),
'R': numpy.array([0.0, 0.5, 0.5]),
'S': numpy.array([0.0, 0.5, 0.0]),
'T': numpy.array([-0.5, 0.5, 0.5]),
'X': numpy.array([zeta, zeta, 0.0]),
'X_1': numpy.array([-zeta, 1 - zeta, 0.0]),
'Y': numpy.array([-0.5, 0.5, 0.0]),
'Z': numpy.array([0.0, 0.0, 0.5])
}
path = [
r"\Gamma", "X", "S", "R", "A", "Z", r"\Gamma", "Y", "X_1", "A_1",
"T", "Y", "Z", "T"
]
return {'kpoints': kpoints, 'path': path}
[docs] def hex(self, is_2d=False):
"""
Generate 'edge' k-points for the hexagonal lattice.
:param bool is_2d: Whether to return 2D path
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' - a list of labels from the 'kpoints'.
"""
self.name = "HEX"
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'A': numpy.array([0.0, 0.0, 0.5]),
'H': numpy.array([1.0 / 3.0, 1.0 / 3.0, 0.5]),
'K': numpy.array([1.0 / 3.0, 1.0 / 3.0, 0.0]),
'L': numpy.array([0.5, 0.0, 0.5]),
'M': numpy.array([0.5, 0.0, 0.0])
}
if is_2d:
self.name += self.TWOD
path = [r"\Gamma", "M", "K", r"\Gamma"]
else:
path = [
r"\Gamma", "M", "K", r"\Gamma", "A", "L", "H", "A", "L", "M",
"K", "H"
]
return {'kpoints': kpoints, 'path': path}
[docs] def rhl1(self, alpha):
"""
Generate 'edge' k-points for the rhombohedral lattice with alpha < 90.
:type alpha: float
:param alpha: Alpha cell angle in radians
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "RHL1"
eta = (1 + 4 * cos(alpha)) / (2 + 4 * cos(alpha))
nu = 3.0 / 4.0 - eta / 2.0
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'B': numpy.array([eta, 0.5, 1.0 - eta]),
'B_1': numpy.array([0.5, 1.0 - eta, eta - 1.0]),
'F': numpy.array([0.5, 0.5, 0.0]),
'L': numpy.array([0.5, 0.0, 0.0]),
'L_1': numpy.array([0.0, 0.0, -0.5]),
'P': numpy.array([eta, nu, nu]),
'P_1': numpy.array([1.0 - nu, 1.0 - nu, 1.0 - eta]),
'P_2': numpy.array([nu, nu, eta - 1.0]),
'Q': numpy.array([1.0 - nu, nu, 0.0]),
'X': numpy.array([nu, 0.0, -nu]),
'Z': numpy.array([0.5, 0.5, 0.5])
}
path = [
r"\Gamma", "L", "B_1", "B", "Z", r"\Gamma", "X", "Q", "F", "P_1",
"Z", "L", "P"
]
return {'kpoints': kpoints, 'path': path}
[docs] def rhl2(self, alpha):
"""
Generate 'edge' k-points for the rhombohedral lattice with alpha > 90.
:type alpha: float
:param alpha: Alpha cell angle in radians
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "RHL2"
eta = 1 / (2 * tan(alpha / 2.0))**2
nu = 0.75 - eta / 2.0
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'F': numpy.array([0.5, -0.5, 0.0]),
'L': numpy.array([0.5, 0.0, 0.0]),
'P': numpy.array([1 - nu, -nu, 1 - nu]),
'P_1': numpy.array([nu, nu - 1.0, nu - 1.0]),
'Q': numpy.array([eta, eta, eta]),
'Q_1': numpy.array([1.0 - eta, -eta, -eta]),
'Z': numpy.array([0.5, -0.5, 0.5])
}
path = [
r"\Gamma", "P", "Z", "Q", r"\Gamma", "F", "P_1", "Q_1", "L", "Z"
]
return {'kpoints': kpoints, 'path': path}
[docs] def mcl(self, b, c, beta):
"""
Generate 'edge' k-points for the monoclinic lattice.
:type b: float
:param b: B length
:type c: float
:param c: C length
:type beta: float
:param beta: Beta cell angle in radians
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "MCL"
eta = (1 - b * cos(beta) / c) / (2 * sin(beta)**2)
nu = 0.5 - eta * c * cos(beta) / b
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'A': numpy.array([0.5, 0.5, 0.0]),
'C': numpy.array([0.0, 0.5, 0.5]),
'D': numpy.array([0.5, 0.0, 0.5]),
'D_1': numpy.array([0.5, 0.5, -0.5]),
'E': numpy.array([0.5, 0.5, 0.5]),
'H': numpy.array([0.0, eta, 1.0 - nu]),
'H_1': numpy.array([0.0, 1.0 - eta, nu]),
'H_2': numpy.array([0.0, eta, -nu]),
'M': numpy.array([0.5, eta, 1.0 - nu]),
'M_1': numpy.array([0.5, 1 - eta, nu]),
'M_2': numpy.array([0.5, 1 - eta, nu]),
'X': numpy.array([0.0, 0.5, 0.0]),
'Y': numpy.array([0.0, 0.0, 0.5]),
'Y_1': numpy.array([0.0, 0.0, -0.5]),
'Z': numpy.array([0.5, 0.0, 0.0])
}
path = [
r"\Gamma", "Y", "H", "C", "E", "M_1", "A", "X", "H_1", "M", "D",
"Z", "Y", "D"
]
return {'kpoints': kpoints, 'path': path}
[docs] def mclc1(self, a, b, c, alpha):
"""
Generate 'edge' k-points for the monoclinic lattice with C setting.
:type a: float
:param a: A length
:type b: float
:param b: B length
:type c: float
:param c: C length
:type alpha: float
:param alpha: Alpha cell angle in radians
:rtype: dict
:return: keys: 'kpoints' is a dict with labels as keys and coordinates
as values, 'path' is a list of labels from the 'kpoints' dict.
"""
self.name = "MCLC1"
zeta = (2 - b * cos(alpha) / c) / (4 * sin(alpha)**2)
eta = 0.5 + 2 * zeta * c * cos(alpha) / b
psi = 0.75 - a**2 / (4 * b**2 * sin(alpha)**2)
phi = psi + (0.75 - psi) * b * cos(alpha) / c
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'N': numpy.array([0.25, 0.25, 0.0]),
'N_1': numpy.array([0.0, -0.5, 0.0]),
'F': numpy.array([1 - zeta, 1 - zeta, 1 - eta]),
'F_1': numpy.array([zeta, zeta, eta]),
'F_2': numpy.array([-zeta, -zeta, 1 - eta]),
'I': numpy.array([phi, 1 - phi, 0.5]),
'I_1': numpy.array([1 - phi, phi - 1, 0.5]),
'L': numpy.array([0.0, 0.5, 0.5]),
'M': numpy.array([0.25, 0.25, 0.5]),
'X': numpy.array([1 - psi, psi - 1, 0.0]),
'X_1': numpy.array([psi, 1 - psi, 0.0]),
'X_2': numpy.array([psi - 1, -psi, 0.0]),
'Y': numpy.array([0.0, 0.5, 0.0]),
'Y_1': numpy.array([0.0, -0.5, 0.0]),
'Z': numpy.array([0.0, 0.0, 0.5])
}
path = [
r"\Gamma", "Y", "F", "L", "I", "I_1", "Z", "F_1", "Y", "X_1", "X",
r"\Gamma", "N", "M", r"\Gamma"
]
return {'kpoints': kpoints, 'path': path}
[docs] def tria(self, is_2d):
self.name = "TRI1a"
kpoints = {
r'\Gamma': numpy.array([0.0, 0.0, 0.0]),
'L': numpy.array([0.5, 0.5, 0.0]),
'M': numpy.array([0.0, 0.5, 0.5]),
'N': numpy.array([0.5, 0.0, 0.5]),
'R': numpy.array([0.5, 0.5, 0.5]),
'X': numpy.array([0.5, 0.0, 0.0]),
'Y': numpy.array([0.0, 0.5, 0.0]),
'Z': numpy.array([0.0, 0.0, 0.5])
}
if is_2d:
self.name += self.TWOD
path = ["X", r"\Gamma", "Y", "L", r"\Gamma"]
else:
path = [
"X", r"\Gamma", "Y", "L", r"\Gamma", "Z", "N", r"\Gamma", "M",
"R", r"\Gamma"
]
return {'kpoints': kpoints, 'path': path}
[docs]class MemoryEstimator(object):
"""
Class to estimate RAM memory for a PW (QE) job.
"""
COMPLEX_SIZE, REAL_SIZE, INT_SIZE = 16.0, 8.0, 4.0
MBYTE = 1024.0 * 1024.0
# Largest FFT dimension that could be factorized
NFFTX = 2049
[docs] def __init__(self, vecs, ecutwfc, ecutrho, nproc, nks, nspin, nbnd, ntyp,
nmix, lscf):
"""
Initialize MemoryEstimator object and set several attributes.
:type vecs: 3 list of 3 floats
:param vecs: Lattice vectors
:type ecutwfc: float
:param ecutwfc: Wavefunction cutoff (Ry)
:type ecutrho: float
:param ecutrho: Density cutoff (Ry)
:type nproc: int
:param nproc: Number of processors
:type nks: int
:param nks: Number of k-points
:type nspin: int
:param nspin: 1 for closed-shell, 2 for spin-polarized
:type nbnd: int
:param nbnd: Number of bands
:type ntyp: int
:param ntyp: Number of atom types
:type nmix: int
:param nmix: Beta mixing
:type lscf: bool
:param lscf: True if the calculation is scf/relax, False if it is nscf
"""
vecs = numpy.array(vecs)
alat_ang = sqrt(numpy.dot(vecs[0], vecs[0]))
at_vecs = vecs / alat_ang
alat = alat_ang * A2B
tpi = 2.0 * pi
fpi = 2.0 * tpi
tpiba = tpi / alat
tpiba2 = tpiba**2
g_fact = 1.0 # gamma_only == FALSE
omega = abs(numpy.dot(numpy.cross(vecs[0], vecs[1]), vecs[2])) * A2B**3
omega_r = tpi**3 / omega
npwx_g = int(fpi / 3.0 * sqrt(ecutwfc)**3 / omega_r / g_fact)
npwx_l = npwx_g / nproc
nktot = nks * nspin + 1 # io_level == 0
npol = 1.0 # noncollin == FALSE
ram = self.COMPLEX_SIZE * nbnd * npol * npwx_l * nktot
# Skip atomic wavefunctions or wannier
# Skip Hubbard wavefunctions
# Skip hybrid functionals
# Skip if defined(__EXX_ACE)
# Skip Nonlocal pseudopotentials V_NL (beta functions), reciprocal space
bg_vecs = numpy.linalg.inv(at_vecs).transpose()
dual = ecutrho / ecutwfc
doublegrid = dual > 4.0
gcutm = dual * ecutwfc / tpiba2
gcutw = ecutwfc / tpiba2
if doublegrid:
gcutms = 4.0 * ecutwfc / tpiba2
else:
gcutms = gcutm
dfftp = self.realSpaceGridInit(at_vecs, bg_vecs, gcutm)
self.dfftp = [self.getGoodFFTDim(x) for x in dfftp]
self.dfftp_nnr = dfftp[0] * dfftp[1] * dfftp[2]
self.ngm = int(fpi / 3.0 * sqrt(ecutrho)**3 / omega_r / g_fact)
if gcutms == gcutm:
self.dffts = list(self.dfftp)
self.ngms = self.ngm
else:
dffts = self.realSpaceGridInit(at_vecs, bg_vecs, gcutms)
self.dffts = [self.getGoodFFTDim(x) for x in dffts]
self.ngms = int(fpi / 3.0 * sqrt(4.0 * ecutwfc)**3 / omega_r)
self.dffts_nnr = self.dffts[0] * self.dffts[1] * self.dffts[2]
# Charge density and potentials - see scf_type in scf_mod
scf_type_size = (self.COMPLEX_SIZE * self.ngm +
self.REAL_SIZE * self.dfftp_nnr) * nspin
# Skip IF ( dft_is_meta() .or. lxdm )
# rho, v, vnew
ram += 3 * scf_type_size
# vltot, vrs, rho_core, rhog_core, psic, strf, kedtau if needed
ram += self.COMPLEX_SIZE * (
self.dfftp_nnr + self.ngm *
(1 + ntyp)) + self.REAL_SIZE * self.dfftp_nnr * (2 + nspin)
# Skip IF ( dft_is_meta() )
if lscf:
# rhoin
ram += scf_type_size
# see mix_type in scf_mod
mix_type_size = self.COMPLEX_SIZE * self.ngm * nspin
# Skip IF ( dft_is_meta() .or. lxdm )
# df, dv (if kept in memory)
ram += mix_type_size * 2 * nmix
# G-vectors: g, gg, mill, nl, nlm, ig_l2g, igtongl
ram += self.REAL_SIZE * self.ngm * 4 + self.INT_SIZE * self.ngm * 7
# double grid: nls, nlsm
ram += self.INT_SIZE * self.ngms * 2
# now scratch space that raises the "high watermark"
ram1 = self.COMPLEX_SIZE * self.ngm * 27 + \
self.REAL_SIZE * self.dffts_nnr
maxram = ram + ram1
maxram = max(
[maxram, self.INT_SIZE * 7 * self.ngm + self.REAL_SIZE * self.ngm])
self.total_ram = maxram / self.MBYTE
[docs] def realSpaceGridInit(self, at_vecs, bg_vecs, gcutm):
"""
Gets minimal 3D real-space FFT grid.
:type at_vecs: 3 x 3 table of floats
:param at_vecs: Lattice vectors in the units of alat
:type bg_vecs: 3 x 3 table of floats
:param bg_vecs: Reciprocal lattice vectors in the units of 1/alat
:type gcutm: float
:param gctum: Radius of the sphere to fit the FFT grid
:rtype: list of 3 floats
:return: FFT grid sizes
"""
# Modules/griddim.f90 :: realspace_grid_init
dfftp = numpy.array(numpy.sqrt(gcutm * numpy.sum(at_vecs**2, axis=1)),
dtype=int)
# See more in QE/Modules/griddim.f90 :: grid_set
mgvecs = numpy.zeros(3)
ranges = [
range(-dfftp[idx], dfftp[idx] + 1) for idx in range(len(dfftp))
]
dfftps = numpy.array(list(product(*ranges)))
gtens = numpy.array(dfftps).dot(bg_vecs)
gsq = numpy.sum(gtens**2, axis=1)
indices = numpy.where(gsq < gcutm)[0]
if len(indices):
maxes = numpy.abs(dfftps[indices]).max(axis=0)
mgvecs = [max(mgvecs[idx], maxes[idx]) for idx in range(len(maxes))]
dfftp = 2 * numpy.asarray(mgvecs) + 1
return dfftp
[docs] def getGoodFFTDim(self, fft_dim):
"""
Get a good FFT dimension.
:type fft_dim: int
:param fft_dim: FFT dimension
:rtype: int
:return: Good FFT dimension, if exists
"""
good_fft_dim = fft_dim
while not self.isGoodDim(good_fft_dim) and good_fft_dim <= self.NFFTX:
good_fft_dim += 1
if good_fft_dim > self.NFFTX:
warnings.warn('Could not find good FFT dimension for %d.' % fft_dim)
return fft_dim
return good_fft_dim
[docs] def isGoodDim(self, fft_dim):
"""
Check if FFT dimension is good.
:type fft_dim: int
:param fft_dim: FFT dimension
:rtype bool
:return: True, if FFT dimension is good, False otherwise
"""
factors = [2.0, 3.0, 5.0, 7.0, 11.0]
# find the factors of the fft dimension
orig_fft_dim = float(fft_dim)
pwr = [0.0] * len(factors)
# Decompose orig_fft_dim into a power product series.
try:
for idx, fac in enumerate(factors):
maxpwr = int(log(orig_fft_dim) / log(fac)) + 1
for jdx in range(maxpwr):
# Dimension is completely decomposed
if orig_fft_dim == 1:
raise StopIteration
if orig_fft_dim % fac == 0:
orig_fft_dim /= fac
pwr[idx] = pwr[idx] + 1
except StopIteration:
pass
fact = orig_fft_dim * numpy.product(
[x**y for x, y in zip(factors, pwr)])
# Check that fact * orig_fft_dim ends up being the same a fft_dim
if fft_dim != fact:
raise warnings.warn('Could not check if FFT dimension is good: '
'%d, %d' % (fft_dim, fact))
# Check if there are other factors larger than 11
if orig_fft_dim != 1:
# FFT dimension contains factors > 11, this is not good
ret = False
else:
# Skipping __ESSL and __LINUX_ESSL
# Check that there are no factors 7 and 11
ret = pwr[-2] == 0 and pwr[-1] == 0
return ret
[docs]class MagSpecies(object):
"""
Class that defines species with starting magnetization.
"""
[docs] def __init__(self, struct=None):
"""
Initialize MagSpecies class.
"""
self.data = dict()
self.species = dict()
if struct:
self.fromStructure(struct)
[docs] def createUniqueElement(self, element, mag_element):
"""
Fill self.data dict. Keys of the self.data are elements. Values are
dicts with magnetization as key and unique element as value. Unique
element is just the atomic symbol plus (if element has more than one
magnetization value) a unique integer.
Example:
{'C': {0.0: 'C'}, 'H': {0.0: 'H', 0.1: 'H1'}}
self.species is a dict where unique elements are keys and elements are
values. Example (based on the example above):
{'C': 'C', 'H': 'H', 'H1': 'H'}
:type element: str
:param element: Element
:type mag: float
:param mag: Starting magnetization
:type hubb_u: float
:param hubb_u: Hubbard U parameter
:rtype: str
:return: Unique element
"""
mag_dict = self.data.get(element)
if mag_dict:
unique_element = mag_dict.get(mag_element)
if unique_element:
return unique_element
else:
unique_element = '%s%d' % (element, len(mag_dict))
self.data[element][mag_element] = unique_element
self.species[unique_element] = element
return unique_element
else:
self.data[element] = {mag_element: element}
self.species[element] = element
return element
[docs] def getMag(self, element, unique_element):
"""
Get magnetization given element and unique element values.
:type element: str
:param element: Element
:type unique_element: str
:param unique_element: Unique element
:rtype: tuple
:return: Starting magnetization and Hubbard U
:raise ValueError: If element, unique_element combination is not found
"""
mag_dict = self.data.get(element)
for key, val in mag_dict.items():
if val == unique_element:
return key
raise ValueError('Could not find magnetization in data (%s), for '
'element (%s) and unique element (%s).' %
(str(self.data), element, unique_element))
[docs] def fromStructure(self, struct):
"""
Set data from the structure.
:param structure.Structure struct: Input structure
"""
self.data = dict()
self.species = dict()
self.elements = []
for atom in struct.atom:
mag_element = get_mag_hubbu(atom)
element = self.createUniqueElement(atom.element, mag_element)
self.elements.append(element)
self.anums = [self.elements.index(x) for x in self.elements]