# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
"""
Low-level objects providing an abstraction for the objects involved in the calculation.
"""
import collections
import abc
import numpy as np
import pymatgen.core.units as units
from pprint import pformat
from monty.design_patterns import singleton
from monty.collections import AttrDict
from enum import Enum
from monty.json import MSONable
from pymatgen.core.units import ArrayWithUnit
from pymatgen.core.lattice import Lattice
from pymatgen.core.structure import Structure
from pymatgen.util.serialization import pmg_serialize
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from monty.json import MontyEncoder, MontyDecoder
[docs]def lattice_from_abivars(cls=None, *args, **kwargs):
"""
Returns a `Lattice` object from a dictionary
with the Abinit variables `acell` and either `rprim` in Bohr or `angdeg`
If acell is not given, the Abinit default is used i.e. [1,1,1] Bohr
Args:
cls: Lattice class to be instantiated. pymatgen.core.lattice.Lattice if `cls` is None
Example:
lattice_from_abivars(acell=3*[10], rprim=np.eye(3))
"""
cls = Lattice if cls is None else cls
kwargs.update(dict(*args))
d = kwargs
rprim = d.get("rprim", None)
angdeg = d.get("angdeg", None)
acell = d["acell"]
if rprim is not None:
if angdeg is not None:
raise ValueError("angdeg and rprimd are mutually exclusive")
rprim = np.reshape(rprim, (3, 3))
rprimd = [float(acell[i]) * rprim[i] for i in range(3)]
# Call pymatgen constructors (note that pymatgen uses Angstrom instead of Bohr).
return cls(ArrayWithUnit(rprimd, "bohr").to("ang"))
elif angdeg is not None:
angdeg = np.reshape(angdeg, 3)
if np.any(angdeg <= 0.):
raise ValueError("Angles must be > 0 but got %s" % str(angdeg))
if angdeg.sum() >= 360.:
raise ValueError("The sum of angdeg must be lower that 360, angdeg %s" % str(angdeg))
# This code follows the implementation in ingeo.F90
# See also http://www.abinit.org/doc/helpfiles/for-v7.8/input_variables/varbas.html#angdeg
tol12 = 1e-12
pi, sin, cos, sqrt = np.pi, np.sin, np.cos, np.sqrt
rprim = np.zeros((3, 3))
if (abs(angdeg[0] - angdeg[1]) < tol12 and abs(angdeg[1] - angdeg[2]) < tol12 and
abs(angdeg[0] - 90.) + abs(angdeg[1] - 90.) + abs(angdeg[2] - 90) > tol12):
# Treat the case of equal angles (except all right angles):
# generates trigonal symmetry wrt third axis
cosang = cos(pi * angdeg[0] / 180.0)
a2 = 2.0 / 3.0 * (1.0 - cosang)
aa = sqrt(a2)
cc = sqrt(1.0 - a2)
rprim[0, 0] = aa
rprim[0, 1] = 0.0
rprim[0, 2] = cc
rprim[1, 0] = -0.5 * aa
rprim[1, 1] = sqrt(3.0) * 0.5 * aa
rprim[1, 2] = cc
rprim[2, 0] = -0.5 * aa
rprim[2, 1] = -sqrt(3.0) * 0.5 * aa
rprim[2, 2] = cc
else:
# Treat all the other cases
rprim[0, 0] = 1.0
rprim[1, 0] = cos(pi * angdeg[2] / 180.)
rprim[1, 1] = sin(pi * angdeg[2] / 180.)
rprim[2, 0] = cos(pi * angdeg[1] / 180.)
rprim[2, 1] = (cos(pi * angdeg[0] / 180.0) - rprim[1, 0] * rprim[2, 0]) / rprim[1, 1]
rprim[2, 2] = sqrt(1.0 - rprim[2, 0] ** 2 - rprim[2, 1] ** 2)
# Call pymatgen constructors (note that pymatgen uses Angstrom instead of Bohr).
rprimd = [float(acell[i]) * rprim[i] for i in range(3)]
return cls(ArrayWithUnit(rprimd, "bohr").to("ang"))
raise ValueError("Don't know how to construct a Lattice from dict:\n%s" % pformat(d))
[docs]def structure_from_abivars(cls=None, *args, **kwargs):
"""
Build a :class:`Structure` object from a dictionary with ABINIT variables.
Args:
cls: Structure class to be instantiated. pymatgen.core.structure.Structure if cls is None
example:
al_structure = structure_from_abivars(
acell=3*[7.5],
rprim=[0.0, 0.5, 0.5,
0.5, 0.0, 0.5,
0.5, 0.5, 0.0],
typat=1,
xred=[0.0, 0.0, 0.0],
ntypat=1,
znucl=13,
)
`xred` can be replaced with `xcart` or `xangst`.
"""
kwargs.update(dict(*args))
d = kwargs
cls = Structure if cls is None else cls
# lattice = Lattice.from_dict(d, fmt="abivars")
lattice = lattice_from_abivars(**d)
coords, coords_are_cartesian = d.get("xred", None), False
if coords is None:
coords = d.get("xcart", None)
if coords is not None:
if "xangst" in d:
raise ValueError("xangst and xcart are mutually exclusive")
coords = ArrayWithUnit(coords, "bohr").to("ang")
else:
coords = d.get("xangst", None)
coords_are_cartesian = True
if coords is None:
raise ValueError("Cannot extract coordinates from:\n %s" % str(d))
coords = np.reshape(coords, (-1, 3))
znucl_type, typat = d["znucl"], d["typat"]
if not isinstance(znucl_type, collections.abc.Iterable):
znucl_type = [znucl_type]
if not isinstance(typat, collections.abc.Iterable):
typat = [typat]
if len(typat) != len(coords):
raise ValueError("len(typat) != len(coords):\ntypat: %s\ncoords: %s" % (typat, coords))
# Note conversion to int and Fortran --> C indexing
typat = np.array(typat, dtype=np.int)
species = [znucl_type[typ - 1] for typ in typat]
return cls(lattice, species, coords, validate_proximity=False,
to_unit_cell=False, coords_are_cartesian=coords_are_cartesian)
[docs]def structure_to_abivars(structure, **kwargs):
"""
Receives a structure and returns a dictionary with the ABINIT variables.
"""
if not structure.is_ordered:
raise ValueError("""\
Received disordered structure with partial occupancies that cannot be converted into an Abinit input
Please use OrderDisorderedStructureTransformation or EnumerateStructureTransformation
to build an appropriate supercell from partial occupancies or alternatively use the Virtual Crystal Approximation.""")
types_of_specie = structure.types_of_specie
natom = structure.num_sites
znucl_type = [specie.number for specie in types_of_specie]
znucl_atoms = structure.atomic_numbers
typat = np.zeros(natom, np.int)
for atm_idx, site in enumerate(structure):
typat[atm_idx] = types_of_specie.index(site.specie) + 1
rprim = ArrayWithUnit(structure.lattice.matrix, "ang").to("bohr")
angdeg = structure.lattice.angles
xred = np.reshape([site.frac_coords for site in structure], (-1, 3))
# Set small values to zero. This usually happens when the CIF file
# does not give structure parameters with enough digits.
rprim = np.where(np.abs(rprim) > 1e-8, rprim, 0.0)
xred = np.where(np.abs(xred) > 1e-8, xred, 0.0)
# Info on atoms.
d = dict(
natom=natom,
ntypat=len(types_of_specie),
typat=typat,
znucl=znucl_type,
xred=xred,
)
# Add info on the lattice.
# Should we use (rprim, acell) or (angdeg, acell) to specify the lattice?
geomode = kwargs.pop("geomode", "rprim")
if geomode == "automatic":
geomode = "rprim"
if structure.lattice.is_hexagonal: # or structure.lattice.is_rhombohedral
geomode = "angdeg"
angdeg = structure.lattice.angles
# Here one could polish a bit the numerical values if they are not exact.
# Note that in pmg the angles are 12, 20, 01 while in Abinit 12, 02, 01
# One should make sure that the orientation is preserved (see Curtarolo's settings)
if geomode == "rprim":
d.update(
acell=3 * [1.0],
rprim=rprim,
)
elif geomode == "angdeg":
d.update(
acell=ArrayWithUnit(structure.lattice.abc, "ang").to("bohr"),
angdeg=angdeg,
)
else:
raise ValueError("Wrong value for geomode: %s" % geomode)
return d
[docs]def contract(s):
"""
>>> assert contract("1 1 1 2 2 3") == "3*1 2*2 1*3"
>>> assert contract("1 1 3 2 3") == "2*1 1*3 1*2 1*3"
"""
if not s:
return s
tokens = s.split()
old = tokens[0]
count = [[1, old]]
for t in tokens[1:]:
if t == old:
count[-1][0] += 1
else:
old = t
count.append([1, t])
return " ".join("%d*%s" % (c, t) for c, t in count)
[docs]class AbivarAble(metaclass=abc.ABCMeta):
"""
An `AbivarAble` object provides a method `to_abivars`
that returns a dictionary with the abinit variables.
"""
[docs] @abc.abstractmethod
def to_abivars(self):
"""Returns a dictionary with the abinit variables."""
# @abc.abstractmethod
# def from_abivars(cls, vars):
# """Build the object from a dictionary with Abinit variables."""
def __str__(self):
return pformat(self.to_abivars(), indent=1, width=80, depth=None)
def __contains__(self, key):
return key in self.to_abivars()
@singleton
class MandatoryVariable:
"""
Singleton used to tag mandatory variables, just because I can use
the cool syntax: variable is MANDATORY!
"""
@singleton
class DefaultVariable:
"""Singleton used to tag variables that will have the default value"""
MANDATORY = MandatoryVariable()
DEFAULT = DefaultVariable()
[docs]class SpinMode(collections.namedtuple('SpinMode', "mode nsppol nspinor nspden"), AbivarAble, MSONable):
"""
Different configurations of the electron density as implemented in abinit:
One can use as_spinmode to construct the object via SpinMode.as_spinmode
(string) where string can assume the values:
- polarized
- unpolarized
- afm (anti-ferromagnetic)
- spinor (non-collinear magnetism)
- spinor_nomag (non-collinear, no magnetism)
"""
[docs] @classmethod
def as_spinmode(cls, obj):
"""Converts obj into a `SpinMode` instance"""
if isinstance(obj, cls):
return obj
else:
# Assume a string with mode
try:
return _mode2spinvars[obj]
except KeyError:
raise KeyError("Wrong value for spin_mode: %s" % str(obj))
[docs] def to_abivars(self):
return {
"nsppol": self.nsppol,
"nspinor": self.nspinor,
"nspden": self.nspden,
}
[docs] @pmg_serialize
def as_dict(self):
return {k: getattr(self, k) for k in self._fields}
[docs] @classmethod
def from_dict(cls, d):
return cls(**{k: d[k] for k in d if k in cls._fields})
# An handy Multiton
_mode2spinvars = {
"unpolarized": SpinMode("unpolarized", 1, 1, 1),
"polarized": SpinMode("polarized", 2, 1, 2),
"afm": SpinMode("afm", 1, 1, 2),
"spinor": SpinMode("spinor", 1, 2, 4),
"spinor_nomag": SpinMode("spinor_nomag", 1, 2, 1),
}
[docs]class Smearing(AbivarAble, MSONable):
"""
Variables defining the smearing technique. The preferred way to instanciate
a `Smearing` object is via the class method Smearing.as_smearing(string)
"""
#: Mapping string_mode --> occopt
_mode2occopt = {
'nosmearing': 1,
'fermi_dirac': 3,
'marzari4': 4,
'marzari5': 5,
'methfessel': 6,
'gaussian': 7}
def __init__(self, occopt, tsmear):
self.occopt = occopt
self.tsmear = tsmear
def __str__(self):
s = "occopt %d # %s Smearing\n" % (self.occopt, self.mode)
if self.tsmear:
s += 'tsmear %s' % self.tsmear
return s
def __eq__(self, other):
return (self.occopt == other.occopt and
np.allclose(self.tsmear, other.tsmear))
def __ne__(self, other):
return not self == other
def __bool__(self):
return self.mode != "nosmearing"
# py2 old version
__nonzero__ = __bool__
[docs] @classmethod
def as_smearing(cls, obj):
"""
Constructs an instance of `Smearing` from obj. Accepts obj in the form:
* Smearing instance
* "name:tsmear" e.g. "gaussian:0.004" (Hartree units)
* "name:tsmear units" e.g. "gaussian:0.1 eV"
* None --> no smearing
"""
if obj is None:
return Smearing.nosmearing()
if isinstance(obj, cls):
return obj
# obj is a string
if obj == "nosmearing":
return cls.nosmearing()
else:
obj, tsmear = obj.split(":")
obj.strip()
occopt = cls._mode2occopt[obj]
try:
tsmear = float(tsmear)
except ValueError:
tsmear, unit = tsmear.split()
tsmear = units.Energy(float(tsmear), unit).to("Ha")
return cls(occopt, tsmear)
@property
def mode(self):
for (mode_str, occopt) in self._mode2occopt.items():
if occopt == self.occopt:
return mode_str
raise AttributeError("Unknown occopt %s" % self.occopt)
[docs] @staticmethod
def nosmearing():
return Smearing(1, 0.0)
[docs] def to_abivars(self):
if self.mode == "nosmearing":
return {"occopt": 1, "tsmear": 0.0}
else:
return {"occopt": self.occopt, "tsmear": self.tsmear, }
[docs] @pmg_serialize
def as_dict(self):
"""json friendly dict representation of Smearing"""
return {"occopt": self.occopt, "tsmear": self.tsmear}
[docs] @staticmethod
def from_dict(d):
return Smearing(d["occopt"], d["tsmear"])
[docs]class ElectronsAlgorithm(dict, AbivarAble, MSONable):
"""Variables controlling the SCF/NSCF algorithm."""
# None indicates that we use abinit defaults.
_DEFAULT = dict(
iprcell=None, iscf=None, diemac=None, diemix=None, diemixmag=None,
dielam=None, diegap=None, dielng=None, diecut=None, nstep=50)
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
for k in self:
if k not in self._DEFAULT:
raise ValueError("%s: No default value has been provided for "
"key %s" % (self.__class__.__name__, k))
[docs] def to_abivars(self):
return self.copy()
[docs] @pmg_serialize
def as_dict(self):
return self.copy()
[docs] @classmethod
def from_dict(cls, d):
d = d.copy()
d.pop("@module", None)
d.pop("@class", None)
return cls(**d)
[docs]class Electrons(AbivarAble, MSONable):
"""The electronic degrees of freedom"""
def __init__(self, spin_mode="polarized", smearing="fermi_dirac:0.1 eV",
algorithm=None, nband=None, fband=None, charge=0.0, comment=None): # occupancies=None,
"""
Constructor for Electrons object.
Args:
comment: String comment for Electrons
charge: Total charge of the system. Default is 0.
"""
super().__init__()
self.comment = comment
self.smearing = Smearing.as_smearing(smearing)
self.spin_mode = SpinMode.as_spinmode(spin_mode)
self.nband = nband
self.fband = fband
self.charge = charge
self.algorithm = algorithm
@property
def nsppol(self):
return self.spin_mode.nsppol
@property
def nspinor(self):
return self.spin_mode.nspinor
@property
def nspden(self):
return self.spin_mode.nspden
[docs] def as_dict(self):
"json friendly dict representation"
d = {}
d["@module"] = self.__class__.__module__
d["@class"] = self.__class__.__name__
d["spin_mode"] = self.spin_mode.as_dict()
d["smearing"] = self.smearing.as_dict()
d["algorithm"] = self.algorithm.as_dict() if self.algorithm else None
d["nband"] = self.nband
d["fband"] = self.fband
d["charge"] = self.charge
d["comment"] = self.comment
return d
[docs] @classmethod
def from_dict(cls, d):
d = d.copy()
d.pop("@module", None)
d.pop("@class", None)
dec = MontyDecoder()
d["spin_mode"] = dec.process_decoded(d["spin_mode"])
d["smearing"] = dec.process_decoded(d["smearing"])
d["algorithm"] = dec.process_decoded(d["algorithm"]) if d["algorithm"] else None
return cls(**d)
[docs] def to_abivars(self):
abivars = self.spin_mode.to_abivars()
abivars.update({
"nband": self.nband,
"fband": self.fband,
"charge": self.charge,
})
if self.smearing:
abivars.update(self.smearing.to_abivars())
if self.algorithm:
abivars.update(self.algorithm)
# abivars["#comment"] = self.comment
return abivars
[docs]class KSamplingModes(Enum):
monkhorst = 1
path = 2
automatic = 3
[docs]class KSampling(AbivarAble, MSONable):
"""
Input variables defining the K-point sampling.
"""
def __init__(self, mode=KSamplingModes.monkhorst, num_kpts=0,
kpts=((1, 1, 1),),
kpt_shifts=(0.5, 0.5, 0.5),
kpts_weights=None, use_symmetries=True, use_time_reversal=True, chksymbreak=None,
comment=None):
"""
Highly flexible constructor for KSampling objects. The flexibility comes
at the cost of usability and in general, it is recommended that you use
the default constructor only if you know exactly what you are doing and
requires the flexibility. For most usage cases, the object be constructed
far more easily using the convenience static constructors:
#. gamma_only
#. gamma_centered
#. monkhorst
#. monkhorst_automatic
#. path
and it is recommended that you use those.
Args:
mode: Mode for generating k-poits. Use one of the KSamplingModes enum types.
num_kpts: Number of kpoints if mode is "automatic"
Number of division for the sampling of the smallest segment if mode is "path".
Not used for the other modes
kpts: Number of divisions. Even when only a single specification is
required, e.g. in the automatic scheme, the kpts should still
be specified as a 2D array. e.g., [[20]] or [[2,2,2]].
kpt_shifts: Shifts for Kpoints.
use_symmetries: False if spatial symmetries should not be used
to reduce the number of independent k-points.
use_time_reversal: False if time-reversal symmetry should not be used
to reduce the number of independent k-points.
kpts_weights: Optional weights for kpoints. For explicit kpoints.
chksymbreak: Abinit input variable: check whether the BZ sampling preserves the symmetry of the crystal.
comment: String comment for Kpoints
.. note::
The default behavior of the constructor is monkhorst.
"""
if isinstance(mode, str):
mode = KSamplingModes[mode]
super().__init__()
self.mode = mode
self.comment = comment
self.num_kpts = num_kpts
self.kpts = kpts
self.kpt_shifts = kpt_shifts
self.kpts_weights = kpts_weights
self.use_symmetries = use_symmetries
self.use_time_reversal = use_time_reversal
self.chksymbreak = chksymbreak
abivars = {}
if mode == KSamplingModes.monkhorst:
assert num_kpts == 0
ngkpt = np.reshape(kpts, 3)
shiftk = np.reshape(kpt_shifts, (-1, 3))
if use_symmetries and use_time_reversal:
kptopt = 1
if not use_symmetries and use_time_reversal:
kptopt = 2
if not use_symmetries and not use_time_reversal:
kptopt = 3
if use_symmetries and not use_time_reversal:
kptopt = 4
abivars.update({
"ngkpt": ngkpt,
"shiftk": shiftk,
"nshiftk": len(shiftk),
"kptopt": kptopt,
"chksymbreak": chksymbreak,
})
elif mode == KSamplingModes.path:
if num_kpts <= 0:
raise ValueError("For Path mode, num_kpts must be specified and >0")
kptbounds = np.reshape(kpts, (-1, 3))
# print("in path with kptbound: %s " % kptbounds)
abivars.update({
"ndivsm": num_kpts,
"kptbounds": kptbounds,
"kptopt": -len(kptbounds) + 1,
})
elif mode == KSamplingModes.automatic:
kpts = np.reshape(kpts, (-1, 3))
if len(kpts) != num_kpts:
raise ValueError("For Automatic mode, num_kpts must be specified.")
abivars.update({
"kptopt": 0,
"kpt": kpts,
"nkpt": num_kpts,
"kptnrm": np.ones(num_kpts),
"wtk": kpts_weights, # for iscf/=-2, wtk.
"chksymbreak": chksymbreak,
})
else:
raise ValueError("Unknown mode %s" % mode)
self.abivars = abivars
# self.abivars["#comment"] = comment
@property
def is_homogeneous(self):
return self.mode not in ["path"]
[docs] @classmethod
def gamma_only(cls):
"""Gamma-only sampling"""
return cls(kpt_shifts=(0.0, 0.0, 0.0), comment="Gamma-only sampling")
[docs] @classmethod
def gamma_centered(cls, kpts=(1, 1, 1), use_symmetries=True, use_time_reversal=True):
"""
Convenient static constructor for an automatic Gamma centered Kpoint grid.
Args:
kpts: Subdivisions N_1, N_2 and N_3 along reciprocal lattice vectors.
use_symmetries: False if spatial symmetries should not be used
to reduce the number of independent k-points.
use_time_reversal: False if time-reversal symmetry should not be used
to reduce the number of independent k-points.
Returns:
:class:`KSampling` object.
"""
return cls(kpts=[kpts], kpt_shifts=(0.0, 0.0, 0.0),
use_symmetries=use_symmetries, use_time_reversal=use_time_reversal,
comment="gamma-centered mode")
[docs] @classmethod
def monkhorst(cls, ngkpt, shiftk=(0.5, 0.5, 0.5), chksymbreak=None, use_symmetries=True,
use_time_reversal=True, comment=None):
"""
Convenient static constructor for a Monkhorst-Pack mesh.
Args:
ngkpt: Subdivisions N_1, N_2 and N_3 along reciprocal lattice vectors.
shiftk: Shift to be applied to the kpoints.
use_symmetries: Use spatial symmetries to reduce the number of k-points.
use_time_reversal: Use time-reversal symmetry to reduce the number of k-points.
Returns:
:class:`KSampling` object.
"""
return cls(
kpts=[ngkpt], kpt_shifts=shiftk,
use_symmetries=use_symmetries, use_time_reversal=use_time_reversal, chksymbreak=chksymbreak,
comment=comment if comment else "Monkhorst-Pack scheme with user-specified shiftk")
[docs] @classmethod
def monkhorst_automatic(cls, structure, ngkpt,
use_symmetries=True, use_time_reversal=True, chksymbreak=None, comment=None):
"""
Convenient static constructor for an automatic Monkhorst-Pack mesh.
Args:
structure: :class:`Structure` object.
ngkpt: Subdivisions N_1, N_2 and N_3 along reciprocal lattice vectors.
use_symmetries: Use spatial symmetries to reduce the number of k-points.
use_time_reversal: Use time-reversal symmetry to reduce the number of k-points.
Returns:
:class:`KSampling` object.
"""
sg = SpacegroupAnalyzer(structure)
# sg.get_crystal_system()
# sg.get_point_group_symbol()
# TODO
nshiftk = 1
# shiftk = 3*(0.5,) # this is the default
shiftk = 3 * (0.5,)
# if lattice.ishexagonal:
# elif lattice.isbcc
# elif lattice.isfcc
return cls.monkhorst(
ngkpt, shiftk=shiftk, use_symmetries=use_symmetries, use_time_reversal=use_time_reversal,
chksymbreak=chksymbreak, comment=comment if comment else "Automatic Monkhorst-Pack scheme")
@classmethod
def _path(cls, ndivsm, structure=None, kpath_bounds=None, comment=None):
"""
Static constructor for path in k-space.
Args:
structure: :class:`Structure` object.
kpath_bounds: List with the reduced coordinates of the k-points defining the path.
ndivsm: Number of division for the smallest segment.
comment: Comment string.
Returns:
:class:`KSampling` object.
"""
if kpath_bounds is None:
# Compute the boundaries from the input structure.
from pymatgen.symmetry.bandstructure import HighSymmKpath
sp = HighSymmKpath(structure)
# Flat the array since "path" is a a list of lists!
kpath_labels = []
for labels in sp.kpath["path"]:
kpath_labels.extend(labels)
kpath_bounds = []
for label in kpath_labels:
red_coord = sp.kpath["kpoints"][label]
# print("label %s, red_coord %s" % (label, red_coord))
kpath_bounds.append(red_coord)
return cls(mode=KSamplingModes.path, num_kpts=ndivsm, kpts=kpath_bounds,
comment=comment if comment else "K-Path scheme")
[docs] @classmethod
def path_from_structure(cls, ndivsm, structure):
"""See _path for the meaning of the variables"""
return cls._path(ndivsm, structure=structure, comment="K-path generated automatically from structure")
[docs] @classmethod
def explicit_path(cls, ndivsm, kpath_bounds):
"""See _path for the meaning of the variables"""
return cls._path(ndivsm, kpath_bounds=kpath_bounds, comment="Explicit K-path")
[docs] @classmethod
def automatic_density(cls, structure, kppa, chksymbreak=None, use_symmetries=True, use_time_reversal=True,
shifts=(0.5, 0.5, 0.5)):
"""
Returns an automatic Kpoint object based on a structure and a kpoint
density. Uses Gamma centered meshes for hexagonal cells and Monkhorst-Pack grids otherwise.
Algorithm:
Uses a simple approach scaling the number of divisions along each
reciprocal lattice vector proportional to its length.
Args:
structure: Input structure
kppa: Grid density
"""
lattice = structure.lattice
lengths = lattice.abc
shifts = np.reshape(shifts, (-1, 3))
ngrid = kppa / structure.num_sites / len(shifts)
mult = (ngrid * lengths[0] * lengths[1] * lengths[2]) ** (1 / 3.)
num_div = [int(round(1.0 / lengths[i] * mult)) for i in range(3)]
# ensure that num_div[i] > 0
num_div = [i if i > 0 else 1 for i in num_div]
angles = lattice.angles
hex_angle_tol = 5 # in degrees
hex_length_tol = 0.01 # in angstroms
right_angles = [i for i in range(3) if abs(angles[i] - 90) < hex_angle_tol]
hex_angles = [i for i in range(3)
if abs(angles[i] - 60) < hex_angle_tol or
abs(angles[i] - 120) < hex_angle_tol]
is_hexagonal = (len(right_angles) == 2 and len(hex_angles) == 1
and abs(lengths[right_angles[0]] -
lengths[right_angles[1]]) < hex_length_tol)
# style = KSamplingModes.gamma
# if not is_hexagonal:
# num_div = [i + i % 2 for i in num_div]
# style = KSamplingModes.monkhorst
comment = "pymatge.io.abinit generated KPOINTS with grid density = " + "{} / atom".format(kppa)
return cls(
mode="monkhorst", num_kpts=0, kpts=[num_div], kpt_shifts=shifts,
use_symmetries=use_symmetries, use_time_reversal=use_time_reversal, chksymbreak=chksymbreak,
comment=comment)
[docs] def to_abivars(self):
return self.abivars
[docs] def as_dict(self):
enc = MontyEncoder()
return {'mode': self.mode.name, 'comment': self.comment,
'num_kpts': self.num_kpts,
'kpts': enc.default(np.array(self.kpts)), 'kpt_shifts': self.kpt_shifts,
'kpts_weights': self.kpts_weights, 'use_symmetries': self.use_symmetries,
'use_time_reversal': self.use_time_reversal, 'chksymbreak': self.chksymbreak,
'@module': self.__class__.__module__, '@class': self.__class__.__name__}
[docs] @classmethod
def from_dict(cls, d):
d = d.copy()
d.pop('@module', None)
d.pop('@class', None)
dec = MontyDecoder()
d['kpts'] = dec.process_decoded(d['kpts'])
return cls(**d)
[docs]class Constraints(AbivarAble):
"""This object defines the constraints for structural relaxation"""
[docs] def to_abivars(self):
raise NotImplementedError("")
[docs]class RelaxationMethod(AbivarAble, MSONable):
"""
This object stores the variables for the (constrained) structural optimization
ionmov and optcell specify the type of relaxation.
The other variables are optional and their use depend on ionmov and optcell.
A None value indicates that we use abinit default. Default values can
be modified by passing them to the constructor.
The set of variables are constructed in to_abivars depending on ionmov and optcell.
"""
_default_vars = {
"ionmov": MANDATORY,
"optcell": MANDATORY,
"ntime": 80,
"dilatmx": 1.05,
"ecutsm": 0.5,
"strfact": None,
"tolmxf": None,
"strtarget": None,
"atoms_constraints": {}, # Constraints are stored in a dictionary. {} means if no constraint is enforced.
}
IONMOV_DEFAULT = 3
OPTCELL_DEFAULT = 2
def __init__(self, *args, **kwargs):
# Initialize abivars with the default values.
self.abivars = self._default_vars
# Overwrite the keys with the args and kwargs passed to constructor.
self.abivars.update(*args, **kwargs)
self.abivars = AttrDict(self.abivars)
for k in self.abivars:
if k not in self._default_vars:
raise ValueError("%s: No default value has been provided for key %s" % (self.__class__.__name__, k))
for k in self.abivars:
if k is MANDATORY:
raise ValueError("%s: No default value has been provided for the mandatory key %s" %
(self.__class__.__name__, k))
[docs] @classmethod
def atoms_only(cls, atoms_constraints=None):
if atoms_constraints is None:
return cls(ionmov=cls.IONMOV_DEFAULT, optcell=0)
else:
return cls(ionmov=cls.IONMOV_DEFAULT, optcell=0, atoms_constraints=atoms_constraints)
[docs] @classmethod
def atoms_and_cell(cls, atoms_constraints=None):
if atoms_constraints is None:
return cls(ionmov=cls.IONMOV_DEFAULT, optcell=cls.OPTCELL_DEFAULT)
else:
return cls(ionmov=cls.IONMOV_DEFAULT, optcell=cls.OPTCELL_DEFAULT, atoms_constraints=atoms_constraints)
@property
def move_atoms(self):
"""True if atoms must be moved."""
return self.abivars.ionmov != 0
@property
def move_cell(self):
"""True if lattice parameters must be optimized."""
return self.abivars.optcell != 0
[docs] def to_abivars(self):
"""Returns a dictionary with the abinit variables"""
# These variables are always present.
out_vars = {
"ionmov": self.abivars.ionmov,
"optcell": self.abivars.optcell,
"ntime": self.abivars.ntime,
}
# Atom relaxation.
if self.move_atoms:
out_vars.update({
"tolmxf": self.abivars.tolmxf,
})
if self.abivars.atoms_constraints:
# Add input variables for constrained relaxation.
raise NotImplementedError("")
out_vars.update(self.abivars.atoms_constraints.to_abivars())
# Cell relaxation.
if self.move_cell:
out_vars.update({
"dilatmx": self.abivars.dilatmx,
"ecutsm": self.abivars.ecutsm,
"strfact": self.abivars.strfact,
"strtarget": self.abivars.strtarget,
})
return out_vars
[docs] def as_dict(self):
d = dict(self._default_vars)
d['@module'] = self.__class__.__module__
d['@class'] = self.__class__.__name__
return d
[docs] @classmethod
def from_dict(cls, d):
d = d.copy()
d.pop('@module', None)
d.pop('@class', None)
return cls(**d)
[docs]class PPModelModes(Enum):
noppmodel = 0
godby = 1
hybersten = 2
linden = 3
farid = 4
[docs]class PPModel(AbivarAble, MSONable):
"""
Parameters defining the plasmon-pole technique.
The common way to instanciate a PPModel object is via the class method PPModel.as_ppmodel(string)
"""
[docs] @classmethod
def as_ppmodel(cls, obj):
"""
Constructs an instance of PPModel from obj.
Accepts obj in the form:
* PPmodel instance
* string. e.g "godby:12.3 eV", "linden".
"""
if isinstance(obj, cls):
return obj
# obj is a string
if ":" not in obj:
mode, plasmon_freq = obj, None
else:
# Extract mode and plasmon_freq
mode, plasmon_freq = obj.split(":")
try:
plasmon_freq = float(plasmon_freq)
except ValueError:
plasmon_freq, unit = plasmon_freq.split()
plasmon_freq = units.Energy(float(plasmon_freq), unit).to("Ha")
return cls(mode=mode, plasmon_freq=plasmon_freq)
def __init__(self, mode="godby", plasmon_freq=None):
if isinstance(mode, str):
mode = PPModelModes[mode]
self.mode = mode
self.plasmon_freq = plasmon_freq
def __eq__(self, other):
if other is None:
return False
else:
if self.mode != other.mode:
return False
if self.plasmon_freq is None:
return other.plasmon_freq is None
else:
return np.allclose(self.plasmon_freq, other.plasmon_freq)
def __ne__(self, other):
return not self == other
def __bool__(self):
return self.mode != PPModelModes.noppmodel
# py2 old version
__nonzero__ = __bool__
def __repr__(self):
return "<%s at %s, mode = %s>" % (self.__class__.__name__, id(self),
str(self.mode))
[docs] def to_abivars(self):
if self:
return {"ppmodel": self.mode.value,
"ppmfrq": self.plasmon_freq}
else:
return {}
[docs] @classmethod
def get_noppmodel(cls):
return cls(mode="noppmodel", plasmon_freq=None)
[docs] def as_dict(self):
return {"mode": self.mode.name, "plasmon_freq": self.plasmon_freq,
"@module": self.__class__.__module__,
"@class": self.__class__.__name__}
[docs] @staticmethod
def from_dict(d):
return PPModel(mode=d["mode"], plasmon_freq=d["plasmon_freq"])
[docs]class ModelDielectricFunction(AbivarAble):
"""Model dielectric function used for BSE calculation"""
def __init__(self, mdf_epsinf):
self.mdf_epsinf = mdf_epsinf
[docs] def to_abivars(self):
return {"mdf_epsinf": self.mdf_epsinf}
##########################################################################################
# WORK IN PROGRESS ######################################
##########################################################################################
[docs]class Screening(AbivarAble):
"""
This object defines the parameters used for the
computation of the screening function.
"""
# Approximations used for W
_WTYPES = {
"RPA": 0,
}
# Self-consistecy modes
_SC_MODES = {
"one_shot": 0,
"energy_only": 1,
"wavefunctions": 2,
}
def __init__(self, ecuteps, nband, w_type="RPA", sc_mode="one_shot",
hilbert=None, ecutwfn=None, inclvkb=2):
"""
Args:
ecuteps: Cutoff energy for the screening (Ha units).
nband Number of bands for the Green's function
w_type: Screening type
sc_mode: Self-consistency mode.
hilbert: Instance of :class:`HilbertTransform` defining the parameters for the Hilber transform method.
ecutwfn: Cutoff energy for the wavefunctions (Default: ecutwfn == ecut).
inclvkb: Option for the treatment of the dipole matrix elements (NC pseudos).
"""
if w_type not in self._WTYPES:
raise ValueError("W_TYPE: %s is not supported" % w_type)
if sc_mode not in self._SC_MODES:
raise ValueError("Self-consistecy mode %s is not supported" % sc_mode)
self.ecuteps = ecuteps
self.nband = nband
self.w_type = w_type
self.sc_mode = sc_mode
self.ecutwfn = ecutwfn
self.inclvkb = inclvkb
if hilbert is not None:
raise NotImplementedError("Hilber transform not coded yet")
self.hilbert = hilbert
# Default values (equivalent to those used in Abinit8)
self.gwpara = 2
self.awtr = 1
self.symchi = 1
self.optdriver = 3
@property
def use_hilbert(self):
return hasattr(self, "hilbert")
# @property
# def gwcalctyp(self):
# "Return the value of the gwcalctyp input variable"
# dig0 = str(self._SIGMA_TYPES[self.type])
# dig1 = str(self._SC_MODES[self.sc_mode]
# return dig1.strip() + dig0.strip()
[docs] def to_abivars(self):
"""Returns a dictionary with the abinit variables"""
abivars = {
"ecuteps": self.ecuteps,
"ecutwfn": self.ecutwfn,
"inclvkb": self.inclvkb,
"gwpara": self.gwpara,
"awtr": self.awtr,
"symchi": self.symchi,
"nband": self.nband,
# "gwcalctyp": self.gwcalctyp,
# "fftgw" : self.fftgw,
"optdriver": self.optdriver,
}
# Variables for the Hilber transform.
if self.use_hilbert:
abivars.update(self.hilbert.to_abivars())
return abivars
[docs]class SelfEnergy(AbivarAble):
"""
This object defines the parameters used for the computation of the self-energy.
"""
_SIGMA_TYPES = {
"gw": 0,
"hartree_fock": 5,
"sex": 6,
"cohsex": 7,
"model_gw_ppm": 8,
"model_gw_cd": 9,
}
_SC_MODES = {
"one_shot": 0,
"energy_only": 1,
"wavefunctions": 2,
}
def __init__(self, se_type, sc_mode, nband, ecutsigx, screening,
gw_qprange=1, ppmodel=None, ecuteps=None, ecutwfn=None, gwpara=2):
"""
Args:
se_type: Type of self-energy (str)
sc_mode: Self-consistency mode.
nband: Number of bands for the Green's function
ecutsigx: Cutoff energy for the exchange part of the self-energy (Ha units).
screening: :class:`Screening` instance.
gw_qprange: Option for the automatic selection of k-points and bands for GW corrections.
See Abinit docs for more detail. The default value makes the code computie the
QP energies for all the point in the IBZ and one band above and one band below the Fermi level.
ppmodel: :class:`PPModel` instance with the parameters used for the plasmon-pole technique.
ecuteps: Cutoff energy for the screening (Ha units).
ecutwfn: Cutoff energy for the wavefunctions (Default: ecutwfn == ecut).
"""
if se_type not in self._SIGMA_TYPES:
raise ValueError("SIGMA_TYPE: %s is not supported" % se_type)
if sc_mode not in self._SC_MODES:
raise ValueError("Self-consistecy mode %s is not supported" % sc_mode)
self.type = se_type
self.sc_mode = sc_mode
self.nband = nband
self.ecutsigx = ecutsigx
self.screening = screening
self.gw_qprange = gw_qprange
self.gwpara = gwpara
if ppmodel is not None:
assert not screening.use_hilbert
self.ppmodel = PPModel.as_ppmodel(ppmodel)
self.ecuteps = ecuteps if ecuteps is not None else screening.ecuteps
self.ecutwfn = ecutwfn
self.optdriver = 4
# band_mode in ["gap", "full"]
# if isinstance(kptgw, str) and kptgw == "all":
# self.kptgw = None
# self.nkptgw = None
# else:
# self.kptgw = np.reshape(kptgw, (-1,3))
# self.nkptgw = len(self.kptgw)
# if bdgw is None:
# raise ValueError("bdgw must be specified")
# if isinstance(bdgw, str):
# # TODO add new variable in Abinit so that we can specify
# # an energy interval around the KS gap.
# homo = float(nele) / 2.0
# #self.bdgw =
# else:
# self.bdgw = np.reshape(bdgw, (-1,2))
# self.freq_int = freq_int
@property
def use_ppmodel(self):
"""True if we are using the plasmon-pole approximation."""
return hasattr(self, "ppmodel")
@property
def gwcalctyp(self):
"""Returns the value of the gwcalctyp input variable."""
dig0 = str(self._SIGMA_TYPES[self.type])
dig1 = str(self._SC_MODES[self.sc_mode])
return dig1.strip() + dig0.strip()
@property
def symsigma(self):
"""1 if symmetries can be used to reduce the number of q-points."""
return 1 if self.sc_mode == "one_shot" else 0
[docs] def to_abivars(self):
"""Returns a dictionary with the abinit variables."""
abivars = dict(
gwcalctyp=self.gwcalctyp,
ecuteps=self.ecuteps,
ecutsigx=self.ecutsigx,
symsigma=self.symsigma,
gw_qprange=self.gw_qprange,
gwpara=self.gwpara,
optdriver=self.optdriver,
nband=self.nband
# "ecutwfn" : self.ecutwfn,
# "kptgw" : self.kptgw,
# "nkptgw" : self.nkptgw,
# "bdgw" : self.bdgw,
)
# FIXME: problem with the spin
# assert len(self.bdgw) == self.nkptgw
# ppmodel variables
if self.use_ppmodel:
abivars.update(self.ppmodel.to_abivars())
return abivars
[docs]class ExcHamiltonian(AbivarAble):
"""This object contains the parameters for the solution of the Bethe-Salpeter equation."""
# Types of excitonic Hamiltonian.
_EXC_TYPES = {
"TDA": 0, # Tamm-Dancoff approximation.
"coupling": 1, # Calculation with coupling.
}
# Algorithms used to compute the macroscopic dielectric function
# and/or the exciton wavefunctions.
_ALGO2VAR = {
"direct_diago": 1,
"haydock": 2,
"cg": 3,
}
# Options specifying the treatment of the Coulomb term.
_COULOMB_MODES = [
"diago",
"full",
"model_df"
]
def __init__(self, bs_loband, nband, mbpt_sciss, coulomb_mode, ecuteps, spin_mode="polarized", mdf_epsinf=None,
exc_type="TDA", algo="haydock", with_lf=True, bs_freq_mesh=None, zcut=None, **kwargs):
"""
Args:
bs_loband: Lowest band index (Fortran convention) used in the e-h basis set.
Can be scalar or array of shape (nsppol,). Must be >= 1 and <= nband
nband: Max band index used in the e-h basis set.
mbpt_sciss: Scissors energy in Hartree.
coulomb_mode: Treatment of the Coulomb term.
ecuteps: Cutoff energy for W in Hartree.
mdf_epsinf: Macroscopic dielectric function :math:`\\epsilon_\\inf` used in
the model dielectric function.
exc_type: Approximation used for the BSE Hamiltonian
with_lf: True if local field effects are included <==> exchange term is included
bs_freq_mesh: Frequency mesh for the macroscopic dielectric function (start, stop, step) in Ha.
zcut: Broadening parameter in Ha.
**kwargs:
Extra keywords
"""
spin_mode = SpinMode.as_spinmode(spin_mode)
# We want an array bs_loband(nsppol).
try:
bs_loband = np.reshape(bs_loband, spin_mode.nsppol)
except ValueError:
bs_loband = np.array(spin_mode.nsppol * [int(bs_loband)])
self.bs_loband = bs_loband
self.nband = nband
self.mbpt_sciss = mbpt_sciss
self.coulomb_mode = coulomb_mode
assert coulomb_mode in self._COULOMB_MODES
self.ecuteps = ecuteps
self.mdf_epsinf = mdf_epsinf
self.exc_type = exc_type
assert exc_type in self._EXC_TYPES
self.algo = algo
assert algo in self._ALGO2VAR
self.with_lf = with_lf
# if bs_freq_mesh is not given, abinit will select its own mesh.
self.bs_freq_mesh = np.array(bs_freq_mesh) if bs_freq_mesh is not None else bs_freq_mesh
self.zcut = zcut
self.optdriver = 99
# Extra options.
self.kwargs = kwargs
# if "chksymbreak" not in self.kwargs:
# self.kwargs["chksymbreak"] = 0
# Consistency check
if any(bs_loband < 0):
raise ValueError("bs_loband <= 0 while it is %s" % bs_loband)
if any(bs_loband >= nband):
raise ValueError("bs_loband (%s) >= nband (%s)" % (bs_loband, nband))
@property
def inclvkb(self):
"""Treatment of the dipole matrix element (NC pseudos, default is 2)"""
return self.kwargs.get("inclvkb", 2)
@property
def use_haydock(self):
"""True if we are using the Haydock iterative technique."""
return self.algo == "haydock"
@property
def use_cg(self):
"""True if we are using the conjugate gradient method."""
return self.algo == "cg"
@property
def use_direct_diago(self):
"""True if we are performing the direct diagonalization of the BSE Hamiltonian."""
return self.algo == "direct_diago"
[docs] def to_abivars(self):
"""Returns a dictionary with the abinit variables."""
abivars = dict(
bs_calctype=1,
bs_loband=self.bs_loband,
# nband=self.nband,
mbpt_sciss=self.mbpt_sciss,
ecuteps=self.ecuteps,
bs_algorithm=self._ALGO2VAR[self.algo],
bs_coulomb_term=21,
mdf_epsinf=self.mdf_epsinf,
bs_exchange_term=1 if self.with_lf else 0,
inclvkb=self.inclvkb,
zcut=self.zcut,
bs_freq_mesh=self.bs_freq_mesh,
bs_coupling=self._EXC_TYPES[self.exc_type],
optdriver=self.optdriver,
)
if self.use_haydock:
# FIXME
abivars.update(
bs_haydock_niter=100, # No. of iterations for Haydock
bs_hayd_term=0, # No terminator
bs_haydock_tol=[0.05, 0], # Stopping criteria
)
elif self.use_direct_diago:
raise NotImplementedError("")
elif self.use_cg:
raise NotImplementedError("")
else:
raise ValueError("Unknown algorithm for EXC: %s" % self.algo)
# Add extra kwargs
abivars.update(self.kwargs)
return abivars