scnc
/
abipy
mirror of https://github.com/abinit/abipy.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/develop' into develop

This commit is contained in:
guido 2017-12-19 23:51:25 +01:00
parent 880b330f1b 5b17a5ef67
commit 90d8c3a596
146 changed files with 14704 additions and 1813 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -1,8 +1,9 @@
# abipy files
abinit_vars.pickle
Profile.prof
abipy/data/runs/flow_*
abipy/examples/flows/flow_*
_debug
__abinb_workdir__
abinb_*.ipynb
*.pickle

6
.travis.yml
View File

@ -9,10 +9,6 @@ matrix:
- os: linux
python: 3.6
env: PYTHON_VERSION=3.6
# - os: osx
# python: 2.7
# language: generic
# env: PYTHON_VERSION=2.7
#dist: precise
sudo: false
@ -100,7 +96,7 @@ branches:
- develop
after_success:
- if [[ "${PYTHON_VERSION}" == "2.7" && "${TRAVIS_OS_NAME}" == "linux" ]]; then coveralls; fi
- if [[ "${PYTHON_VERSION}" == "3.6" && "${TRAVIS_OS_NAME}" == "linux" ]]; then coveralls; fi
notifications:
email:

3
abipy/__init__.py
View File

@ -16,8 +16,6 @@ del sys
#-----------------------------------------------------------------------------
from abipy.core import release
#from abipy.profile import abipy_env
#from abipy.htc import Launcher, MassLauncher, AbinitInput
# Release data
__author__ = ''
@ -27,4 +25,3 @@ del author, email
__license__ = release.license
__version__ = release.version

3
abipy/abilab.py
View File

@ -28,6 +28,7 @@ from abipy.flowtk import Pseudo, PseudoTable, Mrgscr, Mrgddb, Mrggkk, Flow, Task
# g0w0_flow, phonon_flow, phonon_conv_flow, nonlinear_coeff_flow)
from abipy.core.release import __version__, min_abinit_version
from abipy.core.globals import enable_notebook, in_notebook, disable_notebook
from abipy.core import restapi
from abipy.core.structure import (Lattice, Structure, StructureModifier, dataframes_from_structures,
mp_match_structure, mp_search, cod_search)
@ -37,7 +38,7 @@ from abipy.htc.input import AbiInput, LdauParams, LexxParams, input_gen
from abipy.abio.robots import Robot
from abipy.abio.inputs import AbinitInput, MultiDataset, AnaddbInput, OpticInput
from abipy.abio.abivars import AbinitInputFile
from abipy.abio.outputs import AbinitLogFile, AbinitOutputFile, OutNcFile #, CubeFile
from abipy.abio.outputs import AbinitLogFile, AbinitOutputFile, OutNcFile, AboRobot #, CubeFile
#from abipy.tools.plotting import DirTreePlotter
from abipy.tools.printing import print_dataframe
from abipy.tools.notebooks import print_source

18
abipy/abio/factories.py
View File

@ -417,7 +417,7 @@ def ion_ioncell_relax_and_ebands_input(structure, pseudos,
def g0w0_with_ppmodel_inputs(structure, pseudos,
kppa, nscf_nband, ecuteps, ecutsigx,
ecut=None, pawecutdg=None,
ecut=None, pawecutdg=None, shifts=(0.5, 0.5, 0.5),
accuracy="normal", spin_mode="polarized", smearing="fermi_dirac:0.1 eV",
ppmodel="godby", charge=0.0, scf_algorithm=None, inclvkb=2, scr_nband=None,
sigma_nband=None, gw_qprange=1):
@ -447,13 +447,14 @@ def g0w0_with_ppmodel_inputs(structure, pseudos,
See Abinit docs for more detail. The default value makes the code compute the
QP energies for all the point in the IBZ and one band above and one band below the Fermi level.
"""
structure = Structure.as_structure(structure)
multi = MultiDataset(structure, pseudos, ndtset=4)
# Set the cutoff energies.
multi.set_vars(_find_ecut_pawecutdg(ecut, pawecutdg, multi.pseudos))
scf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0)
scf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
scf_electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, algorithm=scf_algorithm,
charge=charge, nband=None, fband=None)
@ -464,7 +465,7 @@ def g0w0_with_ppmodel_inputs(structure, pseudos,
multi[0].set_vars(scf_electrons.to_abivars())
multi[0].set_vars(_stopping_criterion("scf", accuracy))
nscf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0)
nscf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
nscf_electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, algorithm={"iscf": -2},
charge=charge, nband=nscf_nband, fband=None)
@ -482,12 +483,11 @@ def g0w0_with_ppmodel_inputs(structure, pseudos,
multi[2].set_vars(nscf_electrons.to_abivars())
multi[2].set_vars(screening.to_abivars())
multi[2].set_vars(_stopping_criterion("screening", accuracy)) # Dummy
#scr_strategy = ScreeningStrategy(scf_strategy, nscf_strategy, screening)
# Sigma.
if sigma_nband is None: sigma_nband = nscf_nband
self_energy = aobj.SelfEnergy("gw", "one_shot", sigma_nband, ecutsigx, screening,
gw_qprange=gw_qprange, ppmodel=ppmodel)
gw_qprange=gw_qprange, ppmodel=ppmodel)
multi[3].set_vars(nscf_ksampling.to_abivars())
multi[3].set_vars(nscf_electrons.to_abivars())
@ -603,7 +603,7 @@ def g0w0_convergence_inputs(structure, pseudos, kppa, nscf_nband, ecuteps, ecuts
scf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=(0, 0, 0))
nscf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=(0, 0, 0))
else:
# this is the original behaviour before the devellopment of the gwwrapper
# this is the original behaviour before the development of the gwwrapper
scf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0)
nscf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0)
@ -613,7 +613,6 @@ def g0w0_convergence_inputs(structure, pseudos, kppa, nscf_nband, ecuteps, ecuts
charge=charge, nband=max(nscf_nband), fband=None)
multi_scf = MultiDataset(structure, pseudos, ndtset=max(1, len(scf_diffs)))
#print(len(scf_diffs))
multi_scf.set_vars(scf_ksampling.to_abivars())
multi_scf.set_vars(scf_electrons.to_abivars())
@ -660,8 +659,7 @@ def g0w0_convergence_inputs(structure, pseudos, kppa, nscf_nband, ecuteps, ecuts
freqremin=None)
scr_inputs = []
sigma_inputs = []
print("ecuteps", ecuteps, "nscf_nband", nscf_nband)
#print("ecuteps", ecuteps, "nscf_nband", nscf_nband)
for response_model in response_models:
for ecuteps_v in ecuteps:
@ -812,7 +810,7 @@ def scf_phonons_inputs(structure, pseudos, kppa,
# Get the qpoints in the IBZ. Note that here we use a q-mesh with ngkpt=(4,4,4) and shiftk=(0,0,0)
# i.e. the same parameters used for the k-mesh in gs_inp.
qpoints = gs_inp.abiget_ibz(ngkpt=(4,4,4), shiftk=(0,0,0), kptopt=1).points
qpoints = gs_inp.abiget_ibz(ngkpt=(4, 4, 4), shiftk=(0, 0, 0), kptopt=1).points
#print("get_ibz qpoints:", qpoints)
# Build the input files for the q-points in the IBZ.

7
abipy/abio/inputs.py
View File

@ -1091,6 +1091,13 @@ class AbinitInput(six.with_metaclass(abc.ABCMeta, AbstractInput, MSONable, Has_S
"""
return self.remove_vars(_IRDVARS, strict=False)
#def pop_relax_vars(self):
# """
# Remove all the `ird*` variables present in self.
# Return dictionary with the variables that have been removed.
# """
# return scf_input.pop_vars(["ionmov", "optcell", "ntime", "dilatmx"])
@property
def scf_tolvar(self):
"""

202
abipy/abio/outputs.py
View File

@ -11,12 +11,14 @@ from monty.string import is_string
from monty.functools import lazy_property
from monty.termcolor import cprint
from pymatgen.core.units import bohr_to_ang
#from abipy.tools.plotting import set_axlims, add_fig_kwargs, get_ax_fig_plt, get_ax3d_fig_plt
from abipy.core.symmetries import AbinitSpaceGroup
from abipy.core.structure import Structure, dataframes_from_structures
from abipy.core.kpoints import has_timrev_from_kptopt
from abipy.core.mixins import TextFile, AbinitNcFile, NotebookWriter
from abipy.abio.inputs import GEOVARS
from abipy.abio.timer import AbinitTimerParser
from abipy.abio.robots import Robot
from abipy.flowtk import EventsParser, NetcdfReader, GroundStateScfCycle, D2DEScfCycle
@ -47,6 +49,13 @@ class AbinitTextFile(TextFile):
class AbinitLogFile(AbinitTextFile, NotebookWriter):
"""Class representing the Abinit log file."""
def to_string(self, verbose=0):
return str(self.events)
def plot(self, **kwargs):
"""Empty placeholder."""
return None
def write_notebook(self, nbpath=None):
"""
Write an ipython notebook to nbpath. If nbpath is None, a temporay file in the current
@ -140,7 +149,7 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
self.final_vars_global, self.final_vars_dataset = None, None
if self.run_completed:
if self.dryrun_mode:
# footer is not preset. Copy values from header.
# footer is not present. Copy values from header.
self.final_vars_global, self.final_vars_dataset = self.initial_vars_global, self.initial_vars_dataset
else:
self.final_vars_global, self.final_vars_dataset = self._parse_variables("footer")
@ -244,7 +253,7 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
global_kptopt = vars_global.get("kptopt", 1)
structures = []
for i in self.datasets.keys():
for i in self.datasets:
# This code breaks down if there are conflicting GEOVARS in globals and dataset.
d = inigeo.copy()
d.update({k: vars_dataset[i][k] for k in GEOVARS if k in vars_dataset[i]})
@ -498,6 +507,7 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
memory_pre = "P This job should need less than"
magic_exit = "------------- Echo of variables that govern the present computation"
filesizes_pre = "_ WF disk file :"
#verbose = 1
def parse_spgline(line):
"""Parse the line with space group info, return dict."""
@ -517,13 +527,14 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
with open(self.filepath, "rt") as fh:
for line in fh:
line = line.strip()
if verbose: print("inblock:", inblock, " at line:", line)
if line.startswith(magic_exit):
break
if not line or line.startswith("===") or line.startswith("---") or line.startswith("Rough estimation"):
if (not line or line.startswith("===") or line.startswith("---")
or line.startswith("Rough estimation") or line.startswith("PAW method is used")):
continue
if verbose: print("inblock:", inblock, " at line:", line)
if line.startswith("DATASET") or line.startswith("Symmetries :"):
# Get dataset index, parse space group and lattice info, init new dims dict.
inblock = 1
@ -552,6 +563,9 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
assert len(mbpos) == 2
dims["wfk_size_mb"] = float(tokens[mbpos[0]])
dims["denpot_size_mb"] = float(tokens[mbpos[1]])
elif line.startswith("Pmy_natom="):
dims.update(my_natom=int(line.replace("Pmy_natom=", "").strip()))
#print("my_natom", dims["my_natom"])
else:
if line and line[0] == "-": line = line[1:]
tokens = grouper(2, line.replace("=", "").split())
@ -625,7 +639,7 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
other_cycle = other.next_gs_scf_cycle()
if other_cycle is None: break
last = (i == len(others) - 1)
fig = other_cycle.plot(axlist=fig.axes, show=show and last)
fig = other_cycle.plot(ax_list=fig.axes, show=show and last)
if last:
fig.tight_layout()
figures.append(fig)
@ -664,7 +678,7 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
other_cycle = other.next_d2de_scf_cycle()
if other_cycle is None: break
last = (i == len(others) - 1)
fig = other_cycle.plot(axlist=fig.axes, show=show and last)
fig = other_cycle.plot(ax_list=fig.axes, show=show and last)
if last:
fig.tight_layout()
figures.append(fig)
@ -693,48 +707,6 @@ class AbinitOutputFile(AbinitTextFile, NotebookWriter):
return self._write_nb_nbpath(nb, nbpath)
class OutNcFile(AbinitNcFile):
"""
Class representing the _OUT.nc file containing the dataset results
produced at the end of the run. The netcdf variables can be accessed
via instance attribute e.g. `outfile.ecut`. Provides integration with ipython.
"""
def __init__(self, filepath):
super(OutNcFile, self).__init__(filepath)
self.reader = NetcdfReader(filepath)
self._varscache= {k: None for k in self.reader.rootgrp.variables}
def __dir__(self):
"""Ipython integration."""
return sorted(list(self._varscache.keys()))
def __getattribute__(self, name):
try:
return super(OutNcFile, self).__getattribute__(name)
except AttributeError:
# Look in self._varscache
varscache = super(OutNcFile, self).__getattribute__("_varscache")
if name not in varscache:
raise AttributeError("Cannot find attribute %s" % name)
reader = super(OutNcFile, self).__getattribute__("reader")
if varscache[name] is None:
varscache[name] = reader.read_value(name)
return varscache[name]
def close(self):
self.reader.close()
def get_allvars(self):
"""
Read all netcdf variables present in the file.
Return dictionary varname --> value
"""
for k, v in self._varscache.items():
if v is not None: continue
self._varscache[k] = self.reader.read_value(k)
return self._varscache
def validate_output_parser(abitests_dir=None, output_files=None):
"""
Validate/test Abinit output parser.
@ -810,3 +782,135 @@ def validate_output_parser(abitests_dir=None, output_files=None):
cprint("All input files successfully parsed!", "green")
return len(errpaths)
class AboRobot(Robot):
"""
This robot analyzes the results contained in multiple Abinit output files.
Can compare dimensions, SCF cycles, analyze timers.
"""
EXT = "abo"
def get_dims_dataframe(self, index=None):
"""
Build and return pandas DataFrame with the dimensions of the calculation.
Args:
index: Index of the dataframe. Use relative paths of files if None.
"""
rows, my_index = [], []
for i, abo in enumerate(self.abifiles):
try:
dims_dataset, spg_dataset = abo.get_dims_spginfo_dataset()
except Exception as exc:
cprint("Exception while trying to get dimensions from %s\n%s" % (abo.relpath, str(exc)), "yellow")
continue
for dtindex, dims in dims_dataset.items():
dims.update({"dtset": dtindex})
rows.append(dims)
my_index.append(abo.relpath if index is None else index[i])
import pandas as pd
return pd.DataFrame(rows, index=my_index, columns=list(rows[0].keys()))
def get_dataframe(self, with_geo=True, with_dims=True, abspath=False, funcs=None):
"""
Return a pandas DataFrame with the most important results.
and the filenames as index.
Args:
with_geo: True if structure info should be added to the dataframe
with_dims: True if dimensions should be added
abspath: True if paths in index should be absolute. Default: Relative to getcwd().
funcs: Function or list of functions to execute to add more data to the DataFrame.
Each function receives a :class:`GsrFile` object and returns a tuple (key, value)
where key is a string with the name of column and value is the value to be inserted.
"""
rows, row_names = [], []
for label, abo in self:
row_names.append(label)
d = OrderedDict()
if with_dims:
dims_dataset, spg_dataset = abo.get_dims_spginfo_dataset()
if len(dims_dataset) > 1:
cprint("Multiple datasets are not supported. ARGH!", "yellow")
d.update(dims_dataset[1])
# Add info on structure.
if with_geo and abo.run_completed:
d.update(abo.final_structure.get_dict4pandas(with_spglib=True))
# Execute functions
if funcs is not None: d.update(self._exec_funcs(funcs, abo))
rows.append(d)
import pandas as pd
row_names = row_names if not abspath else self._to_relpaths(row_names)
return pd.DataFrame(rows, index=row_names, columns=list(rows[0].keys()))
# TODO
#def gridplot_timer(self)
def write_notebook(self, nbpath=None):
"""
Write a jupyter notebook to nbpath. If nbpath is None, a temporay file in the current
working directory is created. Return path to the notebook.
"""
nbformat, nbv, nb = self.get_nbformat_nbv_nb(title=None)
args = [(l, f.filepath) for l, f in self.items()]
nb.cells.extend([
#nbv.new_markdown_cell("# This is a markdown cell"),
nbv.new_code_cell("robot = abilab.AboRobot(*%s)\nrobot.trim_paths()\nrobot" % str(args)),
nbv.new_code_cell("# robot.get_dims_dataframe()"),
nbv.new_code_cell("robot.get_dataframe()"),
])
# Mixins
nb.cells.extend(self.get_baserobot_code_cells())
return self._write_nb_nbpath(nb, nbpath)
class OutNcFile(AbinitNcFile):
"""
Class representing the _OUT.nc file containing the dataset results
produced at the end of the run. The netcdf variables can be accessed
via instance attribute e.g. `outfile.ecut`. Provides integration with ipython.
"""
def __init__(self, filepath):
super(OutNcFile, self).__init__(filepath)
self.reader = NetcdfReader(filepath)
self._varscache= {k: None for k in self.reader.rootgrp.variables}
def __dir__(self):
"""Ipython integration."""
return sorted(list(self._varscache.keys()))
def __getattribute__(self, name):
try:
return super(OutNcFile, self).__getattribute__(name)
except AttributeError:
# Look in self._varscache
varscache = super(OutNcFile, self).__getattribute__("_varscache")
if name not in varscache:
raise AttributeError("Cannot find attribute %s" % name)
reader = super(OutNcFile, self).__getattribute__("reader")
if varscache[name] is None:
varscache[name] = reader.read_value(name)
return varscache[name]
def close(self):
self.reader.close()
def get_allvars(self):
"""
Read all netcdf variables present in the file.
Return dictionary varname --> value
"""
for k, v in self._varscache.items():
if v is not None: continue
self._varscache[k] = self.reader.read_value(k)
return self._varscache

214
abipy/abio/robots.py
View File

@ -10,11 +10,12 @@ import os
import inspect
from collections import OrderedDict, deque
from functools import wraps
from monty.string import is_string, list_strings
from monty.termcolor import cprint
from abipy.core.mixins import NotebookWriter
from abipy.tools.plotting import plot_xy_with_hue
# TODO: Robot for Abinit output files.
class Robot(NotebookWriter):
"""
@ -30,12 +31,11 @@ class Robot(NotebookWriter):
.. note::
__iter__ returns (label, ncfile)
"""
# TODO
# 2) should or ncfile (not __getitem__ returns ncfiles.__getitem__ !!!
# 3) replace ncfiles with files just to be consistent since we have DdbRobot!
__iter__ returns (label, abifile) so:
for label, abifile in self:
print(label)
"""
# filepaths are relative to `start`. None for asbolute paths. This flag is set in trim_paths
start = None
@ -44,15 +44,18 @@ class Robot(NotebookWriter):
Args:
args is a list of tuples (label, filepath)
"""
self._ncfiles, self._do_close = OrderedDict(), OrderedDict()
self._abifiles, self._do_close = OrderedDict(), OrderedDict()
self._exceptions = deque(maxlen=100)
for label, ncfile in args:
self.add_file(label, ncfile)
#@classmethod
#def from_labels_paths(cls, labels_paths)
# return cls(*labels_paths)
@classmethod
def get_supported_extensions(self):
"""List of strings with extensions supported by Robot subclasses."""
# This is needed to have all subclasses.
from abipy.abilab import Robot
return sorted([cls.EXT for cls in Robot.__subclasses__()])
@classmethod
def class_for_ext(cls, ext):
@ -62,8 +65,8 @@ class Robot(NotebookWriter):
return subcls
raise ValueError("Cannot find Robot subclass associated to extension %s\n" % ext +
"The list of supported extensions is:\n%s" %
[cls.EXT for cls in Robot.__subclasses__()])
"The list of supported extensions (case insensitive) is:\n%s" %
str(cls.get_supported_extensions()))
@classmethod
def from_dir(cls, top, walk=True, abspath=False):
@ -94,7 +97,8 @@ class Robot(NotebookWriter):
items = []
for top in list_strings(dirpaths):
items.extend(cls._open_files_in_dir(top, walk))
if not abspath: new.trim_paths(start=os.path.getcwd())
new = cls(*items)
if not abspath: new.trim_paths(start=os.getcwd())
return new
@classmethod
@ -130,34 +134,42 @@ class Robot(NotebookWriter):
for dirpath, dirnames, filenames in os.walk(top):
filenames = [f for f in filenames if cls.class_handles_filename(f)]
for f in filenames:
ncfile = abiopen(os.path.join(dirpath, f))
if ncfile is not None: items.append((ncfile.filepath, ncfile))
abifile = abiopen(os.path.join(dirpath, f))
if abifile is not None: items.append((abifile.filepath, abifile))
else:
filenames = [f for f in os.listdir(top) if cls.class_handles_filename(f)]
for f in filenames:
ncfile = abiopen(os.path.join(top, f))
if ncfile is not None: items.append((ncfile.filepath, ncfile))
abifile = abiopen(os.path.join(top, f))
if abifile is not None: items.append((abifile.filepath, abifile))
return items
@classmethod
def class_handles_filename(cls, filename):
"""True if robot class handles filename."""
return filename.endswith("_" + cls.EXT + ".nc")
return (filename.endswith("_" + cls.EXT + ".nc") or
filename.endswith("." + cls.EXT)) # This for .abo
@classmethod
def from_files(cls, filenames):
"""Build a Robot from a list of `filenames`."""
def from_files(cls, filenames, labels=None):
"""
Build a Robot from a list of `filenames`.
if labels is None, labels are automatically generated from absolute paths.
"""
filenames = list_strings(filenames)
from abipy.abilab import abiopen
filenames = [f for f in filenames if cls.class_handles_filename(f)]
items = []
for f in filenames:
for i, f in enumerate(filenames):
try:
ncfile = abiopen(f)
abifile = abiopen(f)
except Exception:
ncfile = None
if ncfile is not None: items.append((ncfile.filepath, ncfile))
abifile = None
if abifile is not None:
label = abifile.filepath if labels is None else labels[i]
items.append((label, abifile))
return cls(*items)
@classmethod
@ -182,7 +194,7 @@ class Robot(NotebookWriter):
.. code-block:: python
with GsrRobot.from_flow(flow) as robot:
with abilab.GsrRobot.from_flow(flow) as robot:
print(robot)
# That is equivalent to:
@ -234,6 +246,10 @@ class Robot(NotebookWriter):
"""
if nids and node.node_id not in nids: return
filepath = node.outdir.has_abiext(self.EXT)
if not filepath:
# Look in run.abi directory.
filepath = node.wdir.has_abiext(self.EXT)
if filepath:
try:
label = os.path.relpath(filepath)
@ -241,12 +257,9 @@ class Robot(NotebookWriter):
# current working directory may not be defined!
label = filepath
if task_class is not None:
# Filter by task_class (class or string with class name)
if inspect.isclass(task_class) and not isinstance(node, task_class):
return None
if node.__class__.__name__.lower() != task_class.lower():
return None
# Filter by task_class (class or string with class name)
if task_class is not None and not node.isinstance(task_class):
return None
self.add_file(label, filepath)
@ -262,50 +275,50 @@ class Robot(NotebookWriter):
Number of files found.
"""
count = 0
for filepath, ncfile in self.__class__._open_files_in_dir(top, walk):
for filepath, abifile in self.__class__._open_files_in_dir(top, walk):
count += 1
self.add_file(filepath, ncfile)
self.add_file(filepath, abifile)
return count
def add_file(self, label, ncfile, filter_abifile=None):
def add_file(self, label, abifile, filter_abifile=None):
"""
Add a file to the robot with the given label.
Args:
label: String used to identify the file (must be unique, ax exceptions is
raised if label is already present.
ncfile: Specify the file to be added. Accepts strings (filepath) or abipy file-like objects.
abifile: Specify the file to be added. Accepts strings (filepath) or abipy file-like objects.
filter_abifile: Function that receives an `abifile` object and returns
True if the file should be added to the plotter.
"""
if is_string(ncfile):
if is_string(abifile):
from abipy.abilab import abiopen
ncfile = abiopen(ncfile)
if filter_abifile is not None and not filter_abifile(ncfile):
ncfile.close()
abifile = abiopen(abifile)
if filter_abifile is not None and not filter_abifile(abifile):
abifile.close()
return
# Open file here --> have to close it.
self._do_close[ncfile.filepath] = True
self._do_close[abifile.filepath] = True
if label in self._ncfiles:
if label in self._abifiles:
raise ValueError("label %s is already present!" % label)
self._ncfiles[label] = ncfile
self._abifiles[label] = abifile
#def pop_filepath(self, filepath):
# """
# Remove the file with the given `filepath` and close it.
# """
# if label, ncfile in self._ncfiles.items():
# if ncfile.filepath != filepath: continue
# self._ncfiles.pop(label)
# if self._do_close.pop(ncfile.filepath, False):
# if label, abifile in self._abifiles.items():
# if abifile.filepath != filepath: continue
# self._abifiles.pop(label)
# if self._do_close.pop(abifile.filepath, False):
# try:
# ncfile.close()
# abifile.close()
# except Exception as exc:
# print("Exception while closing: ", ncfile.filepath)
# print("Exception while closing: ", abifile.filepath)
# print(exc)
# #raise
@ -316,11 +329,11 @@ class Robot(NotebookWriter):
return [x for x in list_1 if x in set_2]
#def _get_ointersection_i(self, iattrname):
# if len(self.ncfiles) == 0: return []
# values = list(range(getattr(self.ncfiles[0], iattrname)))
# if len(self.ncfiles) == 1: return values
# for ncfile in self.ncfiles[1:]:
# values = self.ordered_intersection(values, range(getattr(ncfile, iattrname)))
# if len(self.abifiles) == 0: return []
# values = list(range(getattr(self.abifiles[0], iattrname)))
# if len(self.abifiles) == 1: return values
# for abifile in self.abifiles[1:]:
# values = self.ordered_intersection(values, range(getattr(abifile, iattrname)))
# return values
@staticmethod
@ -333,13 +346,13 @@ class Robot(NotebookWriter):
"""
Remove file with the given `label` and close it.
"""
if label in self._ncfiles:
ncfile = self._ncfiles.pop(label)
if self._do_close.pop(ncfile.filepath, False):
if label in self._abifiles:
abifile = self._abifiles.pop(label)
if self._do_close.pop(abifile.filepath, False):
try:
ncfile.close()
abifile.close()
except Exception as exc:
print("Exception while closing: ", ncfile.filepath)
print("Exception while closing: ", abifile.filepath)
print(exc)
#raise
@ -350,15 +363,15 @@ class Robot(NotebookWriter):
if len(new_labels) != len(self):
raise ValueError("Robot has %d files while len(new_labels) = %d" % (len(new_labels), len(self)))
old_labels = list(self._ncfiles.keys())
old_ncfiles = self._ncfiles
self._ncfiles = OrderedDict()
old_labels = list(self._abifiles.keys())
old_abifiles = self._abifiles
self._abifiles = OrderedDict()
new2old = OrderedDict()
for old, new in zip(old_labels, new_labels):
print("old [%s] --> new [%s]" % (old, new))
new2old[new] = old
if not dryrun:
self._ncfiles[new] = old_ncfiles[old]
self._abifiles[new] = old_abifiles[old]
return new2old
@ -368,13 +381,13 @@ class Robot(NotebookWriter):
If start is None, os.getcwd() is used. Set `self.start` attribute, return self.start.
"""
self.start = os.getcwd() if start is None else start
old_paths = list(self._ncfiles.keys())
old_paths = list(self._abifiles.keys())
old_new_paths = [(p, os.path.relpath(os.path.abspath(p), start=self.start)) for p in old_paths]
old_ncfiles = self._ncfiles
self._ncfiles = OrderedDict()
old_abifiles = self._abifiles
self._abifiles = OrderedDict()
for old, new in old_new_paths:
self._ncfiles[new] = old_ncfiles[old]
self._abifiles[new] = old_abifiles[old]
return self.start
@ -384,14 +397,14 @@ class Robot(NotebookWriter):
return self._exceptions
def __len__(self):
return len(self._ncfiles)
return len(self._abifiles)
def __iter__(self):
return iter(self._ncfiles.items())
return iter(self._abifiles.items())
def __getitem__(self, key):
# self[key]
return self.ncfiles.__getitem__(key)
return self._abifiles.__getitem__(key)
def __enter__(self):
return self
@ -401,16 +414,16 @@ class Robot(NotebookWriter):
self.close()
def items(self):
return self._ncfiles.items()
return self._abifiles.items()
@property
def labels(self):
return list(self._ncfiles.keys())
return list(self._abifiles.keys())
def get_label_files_str(self):
"""Return string with [label, filepath]."""
from tabulate import tabulate
return tabulate([(label, ncfile.filepath) for label, ncfile in self], headers=["Label", "Path"])
return tabulate([(label, abifile.relpath) for label, abifile in self], headers=["Label", "Relpath"]) + "\n"
def show_files(self, stream=sys.stdout):
"""Show label --> file path"""
@ -422,9 +435,9 @@ class Robot(NotebookWriter):
def to_string(self, verbose=0):
"""String representation."""
lines = ["%s with %d files in memory:\n" % (self.__class__.__name__, len(self.ncfiles))]
lines = ["%s with %d files in memory:\n" % (self.__class__.__name__, len(self.abifiles))]
app = lines.append
for i, f in enumerate(self.ncfiles):
for i, f in enumerate(self.abifiles):
app(f.to_string(verbose=verbose))
app("\n")
@ -432,12 +445,12 @@ class Robot(NotebookWriter):
def _repr_html_(self):
"""Integration with jupyter notebooks."""
return "<ol>\n{}\n</ol>".format("\n".join("<li>%s</li>" % label for label, ncfile in self))
return "<ol>\n{}\n</ol>".format("\n".join("<li>%s</li>" % label for label, abifile in self))
@property
def ncfiles(self):
def abifiles(self):
"""List of netcdf files."""
return list(self._ncfiles.values())
return list(self._abifiles.values())
def is_sortable(self, aname, raise_exc=False):
"""
@ -445,13 +458,13 @@ class Robot(NotebookWriter):
If raise_exc is True, AttributeError with an explicit message is raised.
"""
try:
obj = getattr(self.ncfiles[0], aname)
obj = getattr(self.abifiles[0], aname)
#return not (callable(obj) or hasattr(obj, "__len__"))
return True
except AttributeError:
if not raise_exc: return False
attrs = []
for key, obj in inspect.getmembers(self.ncfiles[0]):
for key, obj in inspect.getmembers(self.abifiles[0]):
# Ignores anything starting with underscore
if key.startswith('_') or callable(obj) or hasattr(obj, "__len__"): continue
attrs.append(key)
@ -467,23 +480,23 @@ Not all entries are sortable (Please select number-like quantities)""" % (self._
def sortby(self, func_or_string, reverse=False):
"""
Sort files in the robot by `func_or_string`
Return list of (label, ncfile, param) tuples where param is obtained via `func_or_string`.
Return list of (label, abifile, param) tuples where param is obtained via `func_or_string`.
Args:
func_or_string: Either None, string, callable defining the quantity to be used for sorting.
If string, it's assumed that the ncfile has an attribute with the same name and getattr is invoked.
If callable, the output of callable(ncfile) is used.
If string, it's assumed that the abifile has an attribute with the same name and getattr is invoked.
If callable, the output of callable(abifile) is used.
If None, no sorting is performed.
reverse: If set to True, then the list elements are sorted as if each comparison were reversed.
"""
if not func_or_string:
return [(label, ncfile, label) for (label, ncfile) in self]
return [(label, abifile, label) for (label, abifile) in self]
elif callable(func_or_string):
items = [(label, ncfile, func_or_string(ncfile)) for (label, ncfile) in self]
items = [(label, abifile, func_or_string(abifile)) for (label, abifile) in self]
else:
# Assume string and attribute with the same name.
self.is_sortable(func_or_string, raise_exc=True)
items = [(label, ncfile, getattr(ncfile, func_or_string)) for (label, ncfile) in self]
items = [(label, abifile, getattr(abifile, func_or_string)) for (label, abifile) in self]
return sorted(items, key=lambda t: t[2], reverse=reverse)
@ -491,12 +504,12 @@ Not all entries are sortable (Please select number-like quantities)""" % (self._
"""
Close all files that have been opened by the Robot.
"""
for ncfile in self.ncfiles:
if self._do_close.pop(ncfile.filepath, False):
for abifile in self.abifiles:
if self._do_close.pop(abifile.filepath, False):
try:
ncfile.close()
abifile.close()
except:
print("Exception while closing: ", ncfile.filepath)
print("Exception while closing: ", abifile.filepath)
print(exc)
@classmethod
@ -520,8 +533,8 @@ Not all entries are sortable (Please select number-like quantities)""" % (self._
for task in obj.iflat_tasks(nids=nids): #, status=obj.S_OK):
open_method = getattr(task, smeth, None)
if open_method is None: continue
ncfile = open_method()
if ncfile is not None: items.append((task.pos_str, ncfile))
abifile = open_method()
if abifile is not None: items.append((task.pos_str, abifile))
return cls(*items)
else:
@ -531,8 +544,8 @@ Not all entries are sortable (Please select number-like quantities)""" % (self._
#def get_attributes(self, attr_name, obj=None, retdict=False):
# od = OrderedDict()
# for label, ncfile in self:
# obj = ncfile if obj is None else getattr(ncfile, obj)
# for label, abifile in self:
# obj = abifile if obj is None else getattr(abifile, obj)
# od[label] = getattr(obj, attr_name)
# if retdict:
# return od
@ -570,12 +583,12 @@ Not all entries are sortable (Please select number-like quantities)""" % (self._
"""
from abipy.core.structure import dataframes_from_structures
if "index" not in kwargs:
index = list(self._ncfiles.keys())
index = list(self._abifiles.keys())
if not abspath: index = self._to_relpaths(index)
kwargs["index"] = index
ncfiles = self.ncfiles if filter_abifile is not None else list(filter(filter_abifile, self.ncfiles))
return dataframes_from_structures(struct_objects=ncfiles, **kwargs)
abifiles = self.abifiles if filter_abifile is not None else list(filter(filter_abifile, self.abifiles))
return dataframes_from_structures(struct_objects=abifiles, **kwargs)
def get_lattice_dataframe(self, **kwargs):
"""Return pandas DataFrame with lattice parameters."""
@ -587,6 +600,15 @@ Not all entries are sortable (Please select number-like quantities)""" % (self._
dfs = self.get_structure_dataframes(**kwargs)
return dfs.coords
##############################################
# Helper functions to plot pandas dataframes #
##############################################
@staticmethod
@wraps(plot_xy_with_hue)
def plot_xy_with_hue(*args, **kwargs):
return plot_xy_with_hue(*args, **kwargs)
def get_baserobot_code_cells(self, title=None):
"""
Return list of notebook cells with calls to methods provides by the baseclass.

5
abipy/abio/tests/test_inputs.py
View File

@ -6,7 +6,7 @@ import numpy as np
import abipy.data as abidata
from abipy import abilab
from abipy.core.testing import AbipyTest, has_abinit
from abipy.core.testing import AbipyTest
from abipy.abio.inputs import *
from abipy.abio.input_tags import *
@ -443,7 +443,7 @@ class TestAbinitInput(AbipyTest):
#####################
# Non-linear methods
####################
if has_abinit('8.3.2'):
if self.has_abinit(version='8.3.2'):
dte_inputs = gs_inp.make_dte_inputs(phonon_pert=True, skip_permutations=True)
print("dte inputs\n", dte_inputs)
assert len(dte_inputs) == 8
@ -729,6 +729,7 @@ class OpticInputTest(AbipyTest):
)
repr(optic_input); str(optic_input)
assert optic_input.to_string(verbose=2)
# TODO
#assert optic_input._repr_html_()
assert optic_input.vars

16
abipy/abio/tests/test_outputs.py
View File

@ -7,7 +7,7 @@ import abipy.data as abidata
from abipy import abilab
from abipy.core.testing import AbipyTest
from abipy.abio.outputs import AbinitOutputFile, AbinitLogFile
from abipy.abio.outputs import AbinitOutputFile, AbinitLogFile, AboRobot
class AbinitLogFileTest(AbipyTest):
@ -17,6 +17,7 @@ class AbinitLogFileTest(AbipyTest):
log_path = abidata.ref_file("refs/abinit.log")
with AbinitLogFile(log_path) as abilog:
repr(abilog); str(abilog)
assert abilog.to_string(verbose=2)
assert len(abilog.events) == 2
if self.has_nbformat():
abilog.write_notebook(nbpath=self.get_tmpname(text=True))
@ -155,3 +156,16 @@ class AbinitOutputTest(AbipyTest):
assert os.path.exists(abitests_dir)
retcode = validate_output_parser(abitests_dir=abitests_dir)
assert retcode == 0
def test_aborobot(self):
"""Testing AboRobot."""
abo_paths = abidata.ref_files("refs/si_ebands/run.abo", "refs/gs_dfpt.abo")
with AboRobot.from_files(abo_paths) as robot:
repr(robot); str(robot)
assert robot.to_string(verbose=2)
assert robot._repr_html_()
dims = robot.get_dims_dataframe()
df = robot.get_dataframe(with_geo=True)
if self.has_nbformat():
robot.write_notebook(nbpath=self.get_tmpname(text=True))

39
abipy/benchmarks/__init__.py
View File

@ -148,23 +148,38 @@ def bench_main(main):
from abipy.abilab import TaskManager
options.manager = TaskManager.as_manager(options.manager)
#if options.tempdir:
# import tempfile
# options.workdir = tempfile.mkdtemp()
# print("Working in temporary directory", options.workdir)
flow = main(options)
if flow is None: return 0
if options.validate:
# Validate inputs and return
retcode = 0
for task in flow.iflat_tasks():
v = task.input.abivalidate()
if v.retcode != 0: cprint(v, color="red")
retcode += v.retcode
print("input validation retcode: %d" % retcode)
return retcode
if options.abivalidate:
print("Validating input files of the flow...")
isok, errors = flow.abivalidate_inputs()
if not isok:
for e in errors:
if e.retcode == 0: continue
lines = e.log_file.readlines()
i = len(lines) - 50 if len(lines) >= 50 else 0
print("Last 50 line from logfile:")
print("".join(lines[i:]))
raise RuntimeError("flow.abivalidate_input failed. See messages above.")
else:
print("Validation succeeded")
return 0
#if options.remove and os.path.isdir(options.workdir):
# print("Removing old directory:", options.workdir)
# import shutil
# shutil.rmtree(options.workdir)
if options.scheduler:
return flow.make_scheduler().start()
return 0
else:
return flow.build_and_pickle_dump()
return wrapper
@ -198,7 +213,7 @@ def build_bench_main_parser():
parser.add_argument("--min-eff", default=None, type=float, help="Minimum parallel efficiency accepted. Default None.")
parser.add_argument('--paw', default=False, action="store_true", help="Run PAW calculation if available")
parser.add_argument('--validate', default=False, action="store_true", help="Validate input files and return")
parser.add_argument("-a", '--abivalidate', default=False, action="store_true", help="Call Abinit to validate input files and return")
parser.add_argument("-i", '--info', default=False, action="store_true", help="Show benchmark info and exit")
parser.add_argument("-r", "--remove", default=False, action="store_true", help="Remove old flow workdir")

4
abipy/benchmarks/_dfpt_phonon_tio2.py
View File

@ -138,9 +138,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/_string5.py
View File

@ -151,9 +151,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/bse_haydock.py
View File

@ -105,9 +105,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/dfpt_autoparal_phonon.py
View File

@ -122,9 +122,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/dfpt_ddk.py
View File

@ -117,9 +117,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/dfpt_phonon.py
View File

@ -118,9 +118,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/eph_al.py
View File

@ -129,9 +129,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_au108.py
View File

@ -222,9 +222,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_autoparal.py
View File

@ -78,9 +78,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_febio3.py
View File

@ -132,9 +132,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_istwfk.py
View File

@ -78,9 +78,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_kpara.py
View File

@ -69,9 +69,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_mpifft.py
View File

@ -78,9 +78,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_paralkgb.py
View File

@ -92,9 +92,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_pureomp.py
View File

@ -66,9 +66,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_ti256.py
View File

@ -372,9 +372,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_uo2.py
View File

@ -225,9 +225,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gs_useylm.py
View File

@ -71,9 +71,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/gsio_paralkgb.py
View File

@ -80,9 +80,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/moldyn1.py
View File

@ -149,9 +149,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/moldyn5.py
View File

@ -225,9 +225,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/rpa_optic.py
View File

@ -127,9 +127,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/scr_hilbert.py
View File

@ -130,9 +130,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/scr_rpa.py
View File

@ -96,9 +96,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/sigma_cd.py
View File

@ -131,9 +131,7 @@ def main(options):
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

4
abipy/benchmarks/sigma_ppm.py
View File

@ -113,9 +113,7 @@ def main(options):
# print doc string and exit.
print(__doc__)
return
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

48
abipy/core/fields.py
View File

@ -109,7 +109,8 @@ class _Field(Has_Structure):
try:
return self.__class__, float(other)
except:
raise TypeError('object of class %s is not an instance of _Field and cannot be converted to float' % other.__class__)
raise TypeError('object of class %s is not an instance of _Field and cannot be converted to float' %
(other.__class__))
if any([self.nspinor != other.nspinor, self.nsppol != other.nsppol, self.nspden != other.nspden,
self.structure != other.structure, self.mesh != other.mesh]):
@ -258,7 +259,7 @@ class _Field(Has_Structure):
else:
return self.datag.std(axis=axis)
def export(self, filename, visu=None):
def export(self, filename, visu=None, verbose=1):
"""
Export the real space data to file filename.
@ -271,6 +272,7 @@ class _Field(Has_Structure):
visualizer that supports the given file format. If visu is not None, an
instance of visu is returned. See :class:`Visualizer` for the list of
applications and formats supported.
verbose: Verbosity level
Returns:
Instance of :class:`Visualizer`
@ -280,16 +282,22 @@ class _Field(Has_Structure):
tokens = filename.strip().split(".")
ext = tokens[-1]
if verbose:
print("tokens", tokens, "ext", ext)
if not tokens[0]: # filename == ".ext" ==> Create temporary file.
import tempfile
filename = tempfile.mkstemp(suffix="." + ext, text=True)[1]
if not tokens[0]:
# filename == ".ext" ==> Create temporary file.
# dir = os.getcwd() is needed when we invoke the method from a notebook.
# nbworkdir in cwd is needed when we invoke the method from a notebook.
from abipy.core.globals import abinb_mkstemp
_, filename = abinb_mkstemp(suffix="." + ext, text=True)
with open(filename, mode="wt") as fh:
if ext == "xsf":
# xcrysden
xsf.xsf_write_structure(fh, self.structure)
xsf.xsf_write_data(fh, self.structure, self.datar, add_replicas=True)
#elif ext == "POSCAR":
else:
raise NotImplementedError("extension %s is not supported." % ext)
@ -298,24 +306,24 @@ class _Field(Has_Structure):
else:
return visu(filename)
def visualize(self, visu_name):
def visualize(self, appname):
"""
Visualize data with visualizer.
See :class:`Visualizer` for the list of applications and formats supported.
"""
visu = Visualizer.from_name(visu_name)
visu = Visualizer.from_name(appname)
# Try to export data to one of the formats supported by the visualizer
# Use a temporary file (note "." + ext)
for ext in visu.supported_extensions():
ext = "." + ext
try:
return self.export(ext, visu=visu)
return self.export(ext, visu=visu)()
except visu.Error:
pass
else:
raise visu.Error("Don't know how to export data for visualizer %s" % visu_name)
raise visu.Error("Don't know how to export data for visualizer %s" % appname)
def get_interpolator(self):
"""
@ -342,7 +350,7 @@ class _Field(Has_Structure):
# raise NotImplementedError("nspinor != 1 not implenented")
@add_fig_kwargs
def plot_line(self, point1, point2, num=200, cartesian=False, ax=None, **kwargs):
def plot_line(self, point1, point2, num=200, cartesian=False, ax=None, fontsize=12, **kwargs):
"""
Plot (interpolated) density/potential in real space along a line defined
by `point1` and `point2`.
@ -357,6 +365,7 @@ class _Field(Has_Structure):
cartesian: By default, `point1` and `point1` are interpreted as points in fractional
coordinates (if not integers). Use True to pass points in cartesian coordinates.
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: legend and title fontsize.
Return:
`matplotlib` figure
@ -377,12 +386,12 @@ class _Field(Has_Structure):
ax.set_xlabel("Distance in Angstrom from %s at %s" % (r.site1.specie.symbol, cs))
ax.set_ylabel(self.latex_label)
if self.nspden > 1:
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@add_fig_kwargs
def plot_line_neighbors(self, site_index, radius, num=200, max_nn=10, **kwargs):
def plot_line_neighbors(self, site_index, radius, num=200, max_nn=10, fontsize=12, **kwargs):
"""
Plot (interpolated) density/potential in real space along the lines connecting
an atom specified by `site_index` and all neighbors within a sphere of given `radius`.
@ -397,6 +406,7 @@ class _Field(Has_Structure):
radius: Radius of the sphere in Angstrom
num: Number of points sampled along the line.
max_nn: By default, only the first `max_nn` neighbors are showed.
fontsize: legend and title fontsize
Return:
`matplotlib` figure
@ -431,14 +441,14 @@ class _Field(Has_Structure):
ax.plot(r.dist, r.values[ispden],
label=latexlabel_ispden(ispden, self.nspden) if i == 0 else None)
ax.set_title(title)
ax.set_title(title, fontsize=fontsize)
ax.grid(True)
if i == nrows - 1:
ax.set_xlabel("Distance from site_index %s [Angstrom]" % site_index)
ax.set_ylabel(self.latex_label)
if self.nspden > 1:
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@ -569,6 +579,7 @@ class Density(_DensityField):
Initialize the all electron core density of the structure from the pseudopotentials *rhoc* files.
For points close to the atoms, the value at the grid point would be defined as the average on a finer grid
in the neghibourhood of the point.
Args:
valence_density: a Density object representing the valence charge density
structure: the structure for which the total core density will be calculated
@ -579,7 +590,7 @@ class Density(_DensityField):
with respect to nelec. The total density will also be rescaled to fit exactly the given number.
tol: tolerance above which the system will raise an exception if the integrated density doesn't sum up
to the value specified in nelec. Default 0.01 (1% error).
method: different methods to perform the calculation
method: different methods to perform the calculation:
get_sites_in_sphere: based on Structure.get_sites_in_sphere.
mesh3d_dist_gridpoints: based on Mesh3D.dist_gridpoints_in_spheres. Generally faster than
get_sites_in_sphere, but tests can be made for specific cases.
@ -854,30 +865,24 @@ class Density(_DensityField):
# gvecs = self.mesh.gvecs
# gwork = self.mesh.zeros().ravel()
# gnorm = self.structure.gnorm(gvec)
# for idx, gg in enumerate(gvecs):
# #gnorm = self.structure.gnorm(gg)
# gnorm = 1.0 # self.structure.gnorm(gg)
# #gg = np.atleast_2d(gg)
# #mv = np.dot(self.structure.gmet, gg.T)
# #norm2 = 2*np.pi * np.dot(gg, mv)
# #gnorm = np.sqrt(norm2)
# #print gg, gnorm
# if idx != 0:
# gwork[idx] = 4*np.pi/gnorm
# else:
# gwork[idx] = 0.0
# new_shape = self.mesh.ndivs
# gwork = np.reshape(gwork, new_shape)
# #gwork = self.mesh.reshape(gwork)
# # FFT to obtain vh in real space.
# vhg = self.total_rhog * gwork
# vhr = self.mesh.fft_g2r(vhg, fg_ishifted=False)
# return vhr, vhg
def export_to_cube(self, filename, spin='total'):
@ -931,7 +936,6 @@ class Density(_DensityField):
if self.nspinor == 1:
if self.nsppol == 1:
data_dict = {"total": myrhor[0]}
if self.nsppol == 2:
data_dict = {"total": myrhor[0] + myrhor[1], "diff": myrhor[0] - myrhor[1]}

80
abipy/core/globals.py Normal file
View File

@ -0,0 +1,80 @@
"""
Global variables used to initialize AbiPy environment in notebooks.
"""
from __future__ import print_function, division, unicode_literals
from monty.termcolor import cprint
import os
import tempfile
__IN_NOTEBOOK = False
def in_notebook():
"""True if we are running inside a jupyter notebook (and enable_notebook has been called)."""
return __IN_NOTEBOOK
def disable_notebook():
"""Set `in_notebook` flag to False."""
global __IN_NOTEBOOK
__IN_NOTEBOOK = False
def enable_notebook(with_seaborn=True):
"""
Set `in_notebook` flag to True and activate seaborn settings for notebooks if `with_seaborn`.
"""
global __IN_NOTEBOOK
__IN_NOTEBOOK = True
# Use seaborn settings for plots (optional)
if with_seaborn:
import seaborn as sns
sns.set(context='notebook', style='darkgrid', palette='deep',
font='sans-serif', font_scale=1, color_codes=False, rc=None)
def get_abinb_workdir():
"""
Return the absolute path of the scratch directory used to produce
and save temporary files when we are runnning inside a jupyter notebook.
.. note:
Due to web-browser policy, files used in the notebook must be within the cwd.
"""
wdir = os.path.join(os.getcwd(), "__abinb_workdir__")
if not os.path.exists(wdir): os.mkdir(wdir)
return wdir
def abinb_mkstemp(force_abinb_workdir=False, use_relpath=False, **kwargs):
"""
Invoke mkstep with kwargs, return the (fd, name) of the temporary file.
kwargs are passed to `mkstemp` except for `dir` if we are inside a jupyter notebook.
Args:
use_abipy_nbworkdir:
use:relpath: Return relative path (os.path.relpath) if True else absolute (default)
Relative paths are required if we are gonna use the temporary file in
notebooks or in web browers.
In this case, the caller is responsbile for calling the function with the correct flag.
.. example:
_, filename = abinb_mkstep(suffix="." + ext, text=True)
"""
if in_notebook() or force_abinb_workdir:
d = kwargs.pop("dir", None)
if d is not None:
cprint("Files should be created inside abipy_nbworkdir if we are inside jupyter or force_abinb_workdir",
"yellow")
fd, path = tempfile.mkstemp(dir=get_abinb_workdir(), **kwargs)
else:
fd, path = tempfile.mkstemp(**kwargs)
if use_relpath:
path = os.path.relpath(path)
return fd, path

32
abipy/core/kpoints.py
View File

@ -494,7 +494,7 @@ def as_kpoints(obj, lattice, weights=None, names=None):
class Kpoint(SlotPickleMixin):
"""
Class defining one k-point.
Class defining one k-point. This object is immutable and can be used as key in dictionaries
"""
__slots__ = [
@ -570,7 +570,7 @@ class Kpoint(SlotPickleMixin):
"""Set the weight of the k-point."""
self._weight = weight
@property
@lazy_property
def cart_coords(self):
"""Cartesian coordinates of the k-point."""
return self.lattice.get_cartesian_coords(self.frac_coords)
@ -586,7 +586,7 @@ class Kpoint(SlotPickleMixin):
if name is not None and name.startswith("\\"): name = "$" + name + "$"
self._name = name
@property
@lazy_property
def on_border(self):
"""
True if the k-point is on the border of the BZ (lattice translations are taken into account).
@ -599,6 +599,10 @@ class Kpoint(SlotPickleMixin):
return "[%+.3f, %+.3f, %+.3f]" % tuple(self.frac_coords)
def __str__(self):
return self.to_string()
def to_string(self, verbose=0):
"""String representation."""
s = "[%+.3f, %+.3f, %+.3f]" % tuple(self.frac_coords)
if self.name is not None: s += ", name: %s" % self.name
if self._weight is not None: s += ", weight: %.3f" % self.weight
@ -651,7 +655,20 @@ class Kpoint(SlotPickleMixin):
return self.__class__(self.frac_coords.copy(), self.lattice.copy(),
weight=self.weight, name=self.name)
@property
def is_gamma(self, allow_umklapp=False, atol=None):
"""
Return True if this the gamma point.
Args:
allow_umklapp: True if umklapp G-vectors are allowed.
atol: Absolute tolerance for k-point comparison (used only if umklapp).
"""
if not allow_umklapp:
return np.all(self.frac_coords == 0.0)
else:
return issamek(self.frac_coords, [0, 0, 0], atol=atol)
@lazy_property
def norm(self):
"""Norm of the kpoint."""
return np.sqrt(np.dot(self.cart_coords, self.cart_coords))
@ -661,7 +678,6 @@ class Kpoint(SlotPickleMixin):
cls = self.__class__
if self.norm > 1e-12:
return cls(self.frac_coords / self.norm, self.lattice, weight=self.weight)
#except ZeroDivisionError:
else:
return cls.gamma(self.lattice, weight=self.weight)
@ -787,7 +803,11 @@ class KpointList(collections.Sequence):
def to_string(self, func=str, title=None, verbose=0):
"""String representation."""
return "\n".join("%d) %s" % (i, func(kpoint)) for i, kpoint in enumerate(self))
lines = []; app = lines.append
if title is not None: app(marquee(title, mark="="))
lines.extend(["%d) %s" % (i, func(kpoint)) for i, kpoint in enumerate(self)])
return "\n".join(lines)
# Sequence protocol.
def __len__(self):

6
abipy/core/mixins.py
View File

@ -233,14 +233,14 @@ class Has_Structure(object):
"""
return self.structure.export(filepath)
def visualize_structure_with(self, visu_name):
def visualize_structure_with(self, appname):
"""
Visualize the crystalline structure with the specified visualizer.
See :class:`Visualizer` for the list of applications and formats supported.
"""
from abipy.iotools.visualizer import Visualizer
visu = Visualizer.from_name(visu_name)
visu = Visualizer.from_name(appname)
for ext in visu.supported_extensions():
ext = "." + ext
@ -249,7 +249,7 @@ class Has_Structure(object):
except visu.Error:
pass
else:
raise visu.Error("Don't know how to export data for visu_name %s" % visu_name)
raise visu.Error("Don't know how to export data for appname %s" % appname)
@six.add_metaclass(abc.ABCMeta)

39
abipy/core/restapi.py
View File

@ -84,7 +84,7 @@ class PhaseDiagramResults(object):
self.mpids = [e.entry_id for e in entries]
# Create phase diagram.
from pymatgen.phasediagram.maker import PhaseDiagram
from pymatgen.analysis.phase_diagram import PhaseDiagram
self.phasediagram = PhaseDiagram(self.entries)
def plot(self, show_unstable=True, show=True):
@ -100,10 +100,7 @@ class PhaseDiagramResults(object):
Return:
plotter object.
"""
try:
from pymatgen.analysis.phase_diagram import PDPlotter
except ImportError:
from pymatgen.phasediagram.plotter import PDPlotter
from pymatgen.analysis.phase_diagram import PDPlotter
plotter = PDPlotter(self.phasediagram, show_unstable=show_unstable)
if show:
plotter.show()
@ -112,13 +109,10 @@ class PhaseDiagramResults(object):
@lazy_property
def table(self):
"""Pandas dataframe with the most important results."""
# TODO: Use PhaseDiagram but enforce new pymatgen first.
from pymatgen.phasediagram.analyzer import PDAnalyzer
pda = PDAnalyzer(self.phasediagram)
rows = []
for e in self.entries:
d = e.structure.get_dict4pandas(with_spglib=True)
decomp, ehull = pda.get_decomp_and_e_above_hull(e)
decomp, ehull = self.phasediagram.get_decomp_and_e_above_hull(e)
rows.append(OrderedDict([
("Materials ID", e.entry_id),
@ -217,6 +211,7 @@ class DatabaseStructures(NotebookWriter):
"""
Print pandas dataframe, structures using format `fmt`, and data to file `file`.
`fmt` is automaticall set to `cif` if structure is disordered.
Set fmt to None or empty string to disable structure output.
"""
print("\n# Found %s structures in %s database (use `verbose` to get further info)\n"
% (len(self.structures), self.dbname), file=file)
@ -224,15 +219,23 @@ class DatabaseStructures(NotebookWriter):
if self.table is not None: print_dataframe(self.table, file=file)
if verbose and self.data is not None: pprint(self.data, stream=file)
for i, structure in enumerate(self.structures):
print(" ", file=file)
print(marquee("%s input for %s" % (fmt, self.ids[i]), mark="#"), file=file)
print("# " + structure.spget_summary(verbose=verbose).replace("\n", "\n# ") + "\n", file=file)
if not structure.is_ordered:
print(structure.convert(fmt="cif"), file=file)
else:
print(structure.convert(fmt=fmt), file=file)
print(2 * "\n", file=file)
# Print structures
print_structures = not (fmt is None or str(fmt) == "None")
if print_structures:
for i, structure in enumerate(self.structures):
print(" ", file=file)
print(marquee("%s input for %s" % (fmt, self.ids[i]), mark="#"), file=file)
print("# " + structure.spget_summary(verbose=verbose).replace("\n", "\n# ") + "\n", file=file)
if not structure.is_ordered:
print(structure.convert(fmt="cif"), file=file)
else:
print(structure.convert(fmt=fmt), file=file)
print(2 * "\n", file=file)
if len(self.structures) > 10:
# Print info again
print("\n# Found %s structures in %s database (use `verbose` to get further info)\n"
% (len(self.structures), self.dbname), file=file)
def write_notebook(self, nbpath=None, title=None):
"""

24
abipy/core/skw.py
View File

@ -560,7 +560,8 @@ class ElectronInterpolator(object):
self._cached_edos[(kmesh, is_shift)] = edos # .copy()
@add_fig_kwargs
def plot_dos_vs_kmeshes(self, kmeshes, is_shift=None, method="gaussian", step=0.1, width=0.2, ax=None, **kwargs):
def plot_dos_vs_kmeshes(self, kmeshes, is_shift=None, method="gaussian", step=0.1, width=0.2,
fontsize=12, ax=None, **kwargs):
"""
Plot (interpolated) DOSes computed with different meshes.
@ -573,6 +574,7 @@ class ElectronInterpolator(object):
step: Energy step (eV) of the linear mesh.
width: Standard deviation (eV) of the gaussian.
ax: matplotlib `Axes` or None if a new figure should be created.
fontsize: Legend and title fontsize.
Returns:
matplotlib figure.
@ -589,24 +591,26 @@ class ElectronInterpolator(object):
ax.grid(True)
ax.set_xlabel("Energy [eV]")
ax.set_ylabel('DOS [states/eV]')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@add_fig_kwargs
def plot_jdosq0_vs_kmeshes(self, kmeshes, is_shift=None, method="gaussian", step=0.1, width=0.2, ax=None, **kwargs):
def plot_jdosq0_vs_kmeshes(self, kmeshes, is_shift=None, method="gaussian", step=0.1, width=0.2,
ax=None, fontsize=12, **kwargs):
"""
Plot (interpolated) Joint DOSes at q=0 computed with different meshes.
Args:
kmeshes: List of kmeshes. Each item is given by three integers with the number of
divisions along the reciprocal primitive axes.
is_shift: three integers (spglib API). When is_shift is not None, the kmesh is shifted along
is_shift: three integers (spglib API). If None, the kmesh is shifted along
the axis in half of adjacent mesh points irrespective of the mesh numbers. None means unshited mesh.
method: String defining the method for the computation of the DOS.
step: Energy step (eV) of the linear mesh.
width: Standard deviation (eV) of the gaussian.
ax: matplotlib `Axes` or None if a new figure should be created.
fontsize: Legend and title fontsize.
Returns:
matplotlib figure.
@ -622,13 +626,13 @@ class ElectronInterpolator(object):
ax.grid(True)
ax.set_xlabel("Energy [eV]")
ax.set_ylabel('JDOS [states/eV]')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@add_fig_kwargs
def plot_nesting_vs_widths(self, widths, kmesh, e0=None, qvertices_names=None,
line_density=20, is_shift=None, ax=None, **kwargs):
line_density=20, is_shift=None, ax=None, fontsize=12, **kwargs):
"""
Plot (interpolated) nesting factor computed with different broadening.
@ -642,6 +646,7 @@ class ElectronInterpolator(object):
is_shift: three integers (spglib API). When is_shift is not None, the kmesh is shifted along
the axis in half of adjacent mesh points irrespective of the mesh numbers. None means unshited mesh.
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: Legend and title fontsize.
Returns:
matplotlib figure.
@ -658,13 +663,13 @@ class ElectronInterpolator(object):
ax.grid(True)
ax.set_ylabel('Nesting factor')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@add_fig_kwargs
def plot_nesting_vs_kmeshes(self, width, kmeshes, e0=None, qvertices_names=None, line_density=20,
is_shift=None, ax=None, **kwargs):
is_shift=None, ax=None, fontsize=12, **kwargs):
"""
Plot (interpolated) nesting factor computed with different k-meshes.
@ -678,6 +683,7 @@ class ElectronInterpolator(object):
is_shift: three integers (spglib API). When is_shift is not None, the kmesh is shifted along
the axis in half of adjacent mesh points irrespective of the mesh numbers. None means unshited mesh.
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: legend and title fontsize.
Returns:
matplotlib figure.
@ -695,7 +701,7 @@ class ElectronInterpolator(object):
ax.grid(True)
ax.set_ylabel('Nesting factor')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig

148
abipy/core/structure.py
View File

@ -9,6 +9,7 @@ import tempfile
import numpy as np
import pickle
import pymatgen
import pymatgen.core.units as pmg_units
from pprint import pprint, pformat
from warnings import warn
@ -18,7 +19,7 @@ from monty.dev import deprecated
from monty.functools import lazy_property
from monty.string import is_string, marquee
from monty.termcolor import cprint
from pymatgen.core.units import ArrayWithUnit
from pymatgen.core.sites import PeriodicSite
from pymatgen.core.lattice import Lattice
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
@ -350,7 +351,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
"""
cart_coords = np.atleast_2d(cart_coords)
molecule = pymatgen.Molecule([p.symbol for p in pseudos], cart_coords)
l = ArrayWithUnit(acell, "bohr").to("ang")
l = pmg_units.ArrayWithUnit(acell, "bohr").to("ang")
new = molecule.get_boxed_structure(l[0], l[1], l[2])
return cls.as_structure(new)
@ -368,16 +369,18 @@ class Structure(pymatgen.Structure, NotebookWriter):
return cls.boxed_molecule([pseudo], cart_coords, acell=acell)
@classmethod
def bcc(cls, a, species, primitive=True, **kwargs):
def bcc(cls, a, species, primitive=True, units="ang", **kwargs):
"""
Build a primitive or a conventional bcc crystal structure.
Args:
a: Lattice parameter in Angstrom.
a: Lattice parameter (Angstrom if units is not given)
species: Chemical species. See __init__ method of :class:`pymatgen.Structue`
primitive: if True a primitive cell will be produced, otherwise a conventional one
units: Units of input lattice parameters e.g. "bohr", "pm"
kwargs: All keyword arguments accepted by :class:`pymatgen.Structue`
"""
a = pmg_units.Length(a, units).to("ang")
if primitive:
lattice = 0.5 * float(a) * np.array([
-1, 1, 1,
@ -395,16 +398,18 @@ class Structure(pymatgen.Structure, NotebookWriter):
return cls(lattice, species, coords=coords, **kwargs)
@classmethod
def fcc(cls, a, species, primitive=True, **kwargs):
def fcc(cls, a, species, primitive=True, units="ang", **kwargs):
"""
Build a primitive or a conventional fcc crystal structure.
Args:
a: Lattice parameter in Angstrom.
a: Lattice parameter (Angstrom if units is not given)
species: Chemical species. See __init__ method of :class:`pymatgen.Structure`
primitive: if True a primitive cell will be produced, otherwise a conventional one
units: Units of input lattice parameters e.g. "bohr", "pm"
kwargs: All keyword arguments accepted by :class:`pymatgen.Structure`
"""
a = pmg_units.Length(a, units).to("ang")
if primitive:
lattice = 0.5 * float(a) * np.array([
0, 1, 1,
@ -422,12 +427,38 @@ class Structure(pymatgen.Structure, NotebookWriter):
return cls(lattice, species, coords=coords, **kwargs)
@classmethod
def rocksalt(cls, a, species, **kwargs):
def zincblende(cls, a, species, units="ang", **kwargs):
"""
Build a primitive zincblende crystal structure.
Args:
a: Lattice parameter (Angstrom if units is not given)
species: Chemical species. See __init__ method of :class:`pymatgen.Structure`
units: Units of input lattice parameters e.g. "bohr", "pm"
kwargs: All keyword arguments accepted by :class:`pymatgen.Structure`
Example::
Structure.zincblende(a, ["Zn", "S"])
"""
a = pmg_units.Length(a, units).to("ang")
lattice = 0.5 * float(a) * np.array([
0, 1, 1,
1, 0, 1,
1, 1, 0])
frac_coords = np.reshape([0, 0, 0, 0.25, 0.25, 0.25], (2, 3))
return cls(lattice, species, frac_coords, coords_are_cartesian=False, **kwargs)
@classmethod
def rocksalt(cls, a, species, units="ang", **kwargs):
"""
Build a primitive fcc crystal structure.
Args:
a: Lattice parameter in Angstrom.
a: Lattice parameter (Angstrom if units is not given)
units: Units of input lattice parameters e.g. "bohr", "pm"
species: Chemical species. See __init__ method of :class:`pymatgen.Structure`
kwargs: All keyword arguments accepted by :class:`pymatgen.Structure`
@ -436,6 +467,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
Structure.rocksalt(a, ["Na", "Cl"])
"""
a = pmg_units.Length(a, units).to("ang")
lattice = 0.5 * float(a) * np.array([
0, 1, 1,
1, 0, 1,
@ -445,10 +477,17 @@ class Structure(pymatgen.Structure, NotebookWriter):
return cls(lattice, species, frac_coords, coords_are_cartesian=False, **kwargs)
@classmethod
def ABO3(cls, a, species, **kwargs):
def ABO3(cls, a, species, units="ang", **kwargs):
"""
Peroviskite structures.
Args:
a: Lattice parameter (Angstrom if units is not given)
species: Chemical species. See __init__ method of :class:`pymatgen.Structure`
units: Units of input lattice parameters e.g. "bohr", "pm"
kwargs: All keyword arguments accepted by :class:`pymatgen.Structure`
"""
a = pmg_units.Length(a, units).to("ang")
lattice = float(a) * np.eye(3)
frac_coords = np.reshape([
0, 0, 0, # A (2a)
@ -592,7 +631,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
Args:
symprec: Symmetry precision used to refine the structure.
angle_tolerance: Tolerance on anglese
angle_tolerance: Tolerance on angles.
if `symprec` is None and `angle_tolerance` is None, no structure refinement is peformed.
primitive (bool): Whether to convert to a primitive cell following
Setyawan, W., & Curtarolo, S. (2010). doi:10.1016/j.commatsci.2010.05.010
@ -683,28 +722,42 @@ class Structure(pymatgen.Structure, NotebookWriter):
return self.lattice.reciprocal_lattice.matrix
raise ValueError("Wrong value for space: %s " % str(space))
def spget_lattice_type(self, symprec=1e-3, angle_tolerance=5):
"""
Call spglib to get the lattice for the structure, e.g., (triclinic,
orthorhombic, cubic, etc.).This is the same than the
crystal system with the exception of the hexagonal/rhombohedral lattice
Args:
symprec: Symmetry precision for distance
angle_tolerance: Tolerance on angles.
Returns:
(str): Lattice type for structure or None if type cannot be detected.
"""
spgan = SpacegroupAnalyzer(self, symprec=symprec, angle_tolerance=angle_tolerance)
return spgan.get_lattice_type()
def spget_equivalent_atoms(self, symprec=1e-3, angle_tolerance=5, printout=False):
"""
Call spglib to find the inequivalent atoms and build symmetry tables.
Args:
symprec: Symmetry precision for distance
angle_tolerance: Tolerance on anglese
symprec: Symmetry precision for distance.
angle_tolerance: Tolerance on angles.
printout: True to print symmetry tables.
Returns:
namedtuple with the following attributes: irred_pos, eqmap, spgdata
Returns: (irred_pos, eqmap, spgdata)
`namedtuple` with the following attributes:
irred_pos: array giving the position of the i-th irred atom in the structure.
The number of irred atoms is len(irred_pos)
eqmap: Mapping irred atom position --> list with positions of symmetrical atoms
spgdata: spglib dataset with additional data reported by spglib.
irred_pos: array giving the position of the i-th irred atom in the structure.
The number of irred atoms is len(irred_pos)
eqmap: Mapping irred atom position --> list with positions of symmetrical atoms
spgdata: spglib dataset with additional data reported by spglib.
Example::
:Example:
for irr_pos in irred_pos:
eqmap[irr_pos] # List of symmetrical positions associated to the irr_pos atom.
"""
spgan = SpacegroupAnalyzer(self, symprec=symprec, angle_tolerance=angle_tolerance)
spgdata = spgan.get_symmetry_dataset()
@ -738,7 +791,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
Args:
symprec: Symmetry precision for distance
angle_tolerance: Tolerance on anglese
angle_tolerance: Tolerance on angles
verbose: Verbosity level.
"""
spgan = SpacegroupAnalyzer(self, symprec=symprec, angle_tolerance=angle_tolerance)
@ -1030,7 +1083,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
symprec (float): Symmetry precision for structure refinement. If
set to 0, no refinement is done. Otherwise, refinement is
performed using spglib with provided precision.
debye_waller_factors ({element symbol: float}): Allows the
debye_waller_factors `({element symbol: float})`: Allows the
specification of Debye-Waller factors. Note that these
factors are temperature dependent.
two_theta_range ([float of length 2]): Tuple for range of
@ -1045,9 +1098,9 @@ class Structure(pymatgen.Structure, NotebookWriter):
xrd.get_xrd_plot(self, two_theta_range=two_theta_range, annotate_peaks=annotate_peaks, ax=ax)
return fig
def export(self, filename, visu=None):
def export(self, filename, visu=None, verbose=1):
"""
Export the crystalline structure on file filename.
Export the crystalline structure to file `filename`.
Args:
filename: String specifying the file path and the file format.
@ -1056,6 +1109,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
visu: `Visualizer` subclass. By default, this method returns the first available
visualizer that supports the given file format. If visu is not None, an
instance of visu is returned. See :class:`Visualizer` for the list of applications and formats supported.
verbose: Verbosity level
Returns: Instance of :class:`Visualizer`
"""
@ -1064,16 +1118,26 @@ class Structure(pymatgen.Structure, NotebookWriter):
tokens = filename.strip().split(".")
ext = tokens[-1]
#print("tokens", tokens, "ext", ext)
#if ext == "POSCAR":
if not tokens[0]:
# filename == ".ext" ==> Create temporary file.
# dir=os.getcwd() is needed when we invoke the method from a notebook.
_, filename = tempfile.mkstemp(suffix="." + ext, dir=os.getcwd(), text=True)
# nbworkdir in cwd is needed when we invoke the method from a notebook.
from abipy.core.globals import abinb_mkstemp
_, rpath = abinb_mkstemp(force_abinb_workdir=False, use_relpath=False,
suffix="." + ext, text=True)
#if abilab.in_notebook():
# _, filename = tempfile.mkstemp(suffix="." + ext, dir=abilab.get_abipy_nbworkdir(), text=True)
#else:
# _, filename = tempfile.mkstemp(suffix="." + ext, text=True)
if ext == "xsf":
# xcrysden
print("Writing data to:", filename)
s = self.to(fmt="xsf", filename=filename)
if ext.lower() in ("xsf", "poscar", "cif"):
if verbose:
print("Writing data to:", filename, "with fmt:", ext.lower())
s = self.to(fmt=ext)
with open(filename, "wt") as fh:
fh.write(s)
if visu is None:
return Visualizer.from_file(filename)
@ -1109,16 +1173,16 @@ class Structure(pymatgen.Structure, NotebookWriter):
from abipy.display import mvtk
return mvtk.plot_structure(self, figure=figure, show=show, **kwargs)
def visualize(self, visu_name="vesta"):
def visualize(self, appname="vesta"):
"""
Visualize the crystalline structure with visualizer.
See :class:`Visualizer` for the list of applications and formats supported.
"""
if visu_name == "vtk": return self.vtkview()
if visu_name == "mayavi": return self.mayaview()
if appname == "vtk": return self.vtkview()
if appname == "mayavi": return self.mayaview()
# Get the Visualizer subclass from the string.
visu = Visualizer.from_name(visu_name)
visu = Visualizer.from_name(appname)
# Try to export data to one of the formats supported by the visualizer
# Use a temporary file (note "." + ext)
@ -1126,10 +1190,11 @@ class Structure(pymatgen.Structure, NotebookWriter):
ext = "." + ext
try:
return self.export(ext, visu=visu)()
except visu.Error:
except visu.Error as exc:
print(exc)
pass
else:
raise visu.Error("Don't know how to export data for %s" % visu_name)
raise visu.Error("Don't know how to export data for %s" % appname)
def convert(self, fmt="cif", **kwargs):
"""
@ -1670,7 +1735,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
pseudo = table.pseudo_with_symbol(site.species_string)
nval += pseudo.Z_val
return nval
return int(nval) if int(nval) == nval else nval
def valence_electrons_per_atom(self, pseudos):
"""
@ -1732,13 +1797,16 @@ def dataframes_from_structures(struct_objects, index=None, with_spglib=True, car
cart_coords: True if the `coords` dataframe should contain Cartesian cordinates
instead of Reduced coordinates.
Return:
Return: (lattice, coords)
namedtuple with two pandas :class:`DataFrame`:
`lattice` contains the lattice parameters,
`coords` the atomic positions.
The list of structures is available in the `structures` entry.
Example::
:Example:
dfs = dataframes_from_structures(files)
dfs.lattice

53
abipy/core/testing.py
View File

@ -15,7 +15,9 @@ import tempfile
import shutil
import unittest
import numpy.testing.utils as nptu
import abipy.data as abidata
from functools import wraps
from monty.os.path import which
from monty.string import is_string
from pymatgen.util.testing import PymatgenTest
@ -146,7 +148,6 @@ def json_read_abinit_input_from_path(json_path):
Read a json file from the absolute path `json_path`, return AbinitInput instance.
"""
from abipy.abio.inputs import AbinitInput
import abipy.data as abidata
with open(json_path, "rt") as fh:
d = json.load(fh)
@ -167,7 +168,6 @@ def input_equality_check(ref_file, input2, rtol=1e-05, atol=1e-08, equal_nan=Fal
we check if all vars are uniquely present in both inputs and if the values are equal (integers, strings)
or almost equal (floats)
"""
def check_int(i, j):
return i != j
@ -246,6 +246,35 @@ def input_equality_check(ref_file, input2, rtol=1e-05, atol=1e-08, equal_nan=Fal
raise AssertionError(msg)
def get_gsinput_si(usepaw=0, as_task=False):
# Build GS input file.
pseudos = abidata.pseudos("14si.pspnc") if usepaw == 0 else data.pseudos("Si.GGA_PBE-JTH-paw.xml")
#silicon = abilab.Structure.zincblende(5.431, ["Si", "Si"], units="ang")
silicon = abidata.cif_file("si.cif")
from abipy.abio.inputs import AbinitInput
scf_input = AbinitInput(silicon, pseudos)
ecut = 6
scf_input.set_vars(
ecut=ecut,
pawecutdg=40,
nband=6,
paral_kgb=0,
iomode=3,
toldfe=1e-9,
)
if usepaw:
scf_input.set_vars(pawecutdg=4 * ecut)
# K-point sampling (shifted)
scf_input.set_autokmesh(nksmall=4)
if not as_task:
return scf_input
else:
from abipy.flowtk.tasks import ScfTask
return ScfTask(scf_input)
class AbipyTest(PymatgenTest):
"""Extends PymatgenTest with Abinit-specific methods """
@ -301,7 +330,6 @@ class AbipyTest(PymatgenTest):
@staticmethod
def get_abistructure_from_abiref(basename):
"""Return an Abipy structure from the basename of one of the reference files."""
import abipy.data as abidata
from abipy.core.structure import Structure
return Structure.as_structure(abidata.ref_file(basename))
@ -400,19 +428,6 @@ class AbipyTest(PymatgenTest):
msg = "This test requires phonopy version %s %s" % (op, version)
raise unittest.SkipTest(msg)
@staticmethod
def skip_if_not_abinit(version=None, op=">="):
"""
Raise SkipTest if the specified version of abinit is not installed.
Use `version` and `op` to ask for a specific version
"""
if not has_abinit(version=version, op=op):
if version is None:
msg = "This test requires abinit"
else:
msg = "This test requires abinit version %s %s" % (op, version)
raise unittest.SkipTest(msg)
@staticmethod
def skip_if_not_pseudodojo():
"""
@ -470,6 +485,7 @@ class AbipyTest(PymatgenTest):
assert not errors
@staticmethod
def abivalidate_flow(self, flow):
"""
Invoke Abinit to test validity of the inputs of a flow
@ -486,6 +502,11 @@ class AbipyTest(PymatgenTest):
print("".join(lines[i:]))
raise RuntimeError("flow.abivalidate_input failed. See messages above.")
@staticmethod
@wraps(get_gsinput_si)
def get_gsinput_si(*args, **kwargs):
return get_gsinput_si(*args, **kwargs)
class AbipyFileTest(AbipyTest):
"""

4
abipy/core/tests/test_fields.py
View File

@ -41,8 +41,8 @@ class TestScalarField(AbipyTest):
self.assert_equal(field.std(space="g"), field.datag.std(axis=0))
field.export(self.get_tmpname(text=True, suffix=".xsf"))
visu = field.visualize("vesta")
assert callable(visu)
#visu = field.visualize(appname="vesta")
#assert callable(visu)
# Field "algebra"
#assert field.datar_xyz.ndim == 4

6
abipy/core/tests/test_kpoints.py
View File

@ -86,6 +86,12 @@ class TestKpoint(AbipyTest):
X = Kpoint([0.5, 0, 0], lattice)
K = Kpoint([1/3, 1/3, 1/3], lattice)
repr(X); str(X)
assert X.to_string(verbose=2)
assert gamma.is_gamma()
assert not pgamma.is_gamma()
assert pgamma.is_gamma(allow_umklapp=True)
assert not X.is_gamma()
# TODO
#assert np.all(np.array(X) == X.frac_coords)

13
abipy/core/tests/test_structure.py
View File

@ -35,7 +35,7 @@ class TestStructure(AbipyTest):
# Export data in Xcrysden format.
#structure.export(self.get_tmpname(text=True, suffix=".xsf"))
#visu = structure.visualize(visu_name="vesta")
#visu = structure.visualize(appname="vesta")
#assert callable(visu)
if self.has_ase():
@ -47,6 +47,7 @@ class TestStructure(AbipyTest):
si = Structure.as_structure(abidata.cif_file("si.cif"))
assert si.formula == "Si2"
assert si.abi_spacegroup is None and not si.has_abi_spacegroup
assert "ntypat" in si.to(fmt="abivars")
spgroup = si.spgset_abi_spacegroup(has_timerev=True)
assert spgroup is not None
@ -113,6 +114,9 @@ class TestStructure(AbipyTest):
assert Specie("Zn", 2) in oxi_znse.composition.elements
assert Specie("Se", -2) in oxi_znse.composition.elements
system = si.spget_lattice_type()
assert system == "cubic"
e = si.spget_equivalent_atoms(printout=True)
assert len(e.irred_pos) == 1
self.assert_equal(e.eqmap[0], [0, 1])
@ -149,6 +153,7 @@ class TestStructure(AbipyTest):
mgb2_cod = Structure.from_cod_id(1526507, primitive=True)
assert mgb2_cod.formula == "Mg1 B2"
assert mgb2_cod.spget_lattice_type() == "hexagonal"
mgb2 = abidata.structure_from_ucell("MgB2")
if self.has_ase():
@ -167,7 +172,8 @@ class TestStructure(AbipyTest):
self.serialize_with_pickle(mgb2)
pseudos = abidata.pseudos("12mg.pspnc", "5b.pspnc")
assert mgb2.num_valence_electrons(pseudos) == 8
nv = mgb2.num_valence_electrons(pseudos)
assert nv == 8 and isinstance(nv , int)
assert mgb2.valence_electrons_per_atom(pseudos) == [2, 3, 3]
self.assert_equal(mgb2.calc_shiftk() , [[0.0, 0.0, 0.5]])
@ -184,6 +190,7 @@ class TestStructure(AbipyTest):
# Function to compute cubic a0 from primitive v0 (depends on struct_type)
vol2a = {"fcc": lambda vol: (4 * vol) ** (1/3.),
"bcc": lambda vol: (2 * vol) ** (1/3.),
"zincblende": lambda vol: (4 * vol) ** (1/3.),
"rocksalt": lambda vol: (4 * vol) ** (1/3.),
"ABO3": lambda vol: vol ** (1/3.),
"hH": lambda vol: (4 * vol) ** (1/3.),
@ -202,6 +209,8 @@ class TestStructure(AbipyTest):
fcc_conv = Structure.fcc(a, ["Si"], primitive=False)
assert len(fcc_conv) == 4
self.assert_almost_equal(a**3, fcc_conv.volume)
zns = Structure.zincblende(a / bohr_to_ang, ["Zn", "S"], units="bohr")
self.assert_almost_equal(a, vol2a["zincblende"](zns.volume))
rock = Structure.rocksalt(a, ["Na", "Cl"])
assert len(rock) == 2
self.assert_almost_equal(a, vol2a["rocksalt"](rock.volume))

2
abipy/data/managers/nic4_manager.yml
View File

@ -6,8 +6,8 @@ hardware: &hardware
mem_per_node: 64Gb
job: &job
#mpi_runner: "mpirun --bind-to none "
mpi_runner: "mpirun"
mpi_runner_options: "--bind-to none"
shell_env:
PATH: "/home/users/g/m/gmatteo/git_repos/abinit/_build_nic4-intel-openmpi-mkl-hdf5.ac/src/98_main:$PATH"
pre_run:

4875
abipy/data/pseudos/As_r.psp8 Normal file
View File

File diff suppressed because it is too large Load Diff

4875
abipy/data/pseudos/Ga-low_r.psp8 Normal file
View File

File diff suppressed because it is too large Load Diff

7
abipy/data/refs/al_eph/_runflow.py
View File

@ -98,11 +98,10 @@ def build_flow(options):
return flow
@abilab.flow_main
@flowtk.flow_main
def main(options):
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":
retcode = main()

10
abipy/data/refs/diamond_sigeph/_runflow.py
View File

@ -113,17 +113,15 @@ def build_flow(options):
flow.allocate()
# EPH does not support autoparal (yet)
for eph_task in eph_work:
eph_task.with_fixed_mpi_omp(1, 1)
#for eph_task in eph_work:
# eph_task.with_fixed_mpi_omp(1, 1)
return flow
@abilab.flow_main
@flowtk.flow_main
def main(options):
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":
retcode = main()

6
abipy/data/refs/gaas_optic/_runflow.py
View File

@ -83,11 +83,9 @@ def build_flow(options, paral_kgb=0):
return flow
@abilab.flow_main
@flowtk.flow_main
def main(options):
flow = build_flow(options)
flow.build_and_pickle_dump()
return flow
return build_flow(options)
if __name__ == "__main__":

278
abipy/dfpt/ddb.py
View File

@ -5,9 +5,11 @@ from __future__ import print_function, division, unicode_literals, absolute_impo
import sys
import os
import tempfile
import itertools
import numpy as np
from collections import OrderedDict
import pandas as pd
from collections import OrderedDict
from six.moves import map, zip, StringIO
from monty.string import marquee
from monty.collections import AttrDict, dict2namedtuple, tree
@ -18,7 +20,7 @@ from abipy.flowtk import NetcdfReader, AnaddbTask
from abipy.core.mixins import TextFile, Has_Structure, NotebookWriter
from abipy.core.symmetries import AbinitSpaceGroup
from abipy.core.structure import Structure
from abipy.core.kpoints import KpointList
from abipy.core.kpoints import KpointList, Kpoint
from abipy.core.tensor import Tensor
from abipy.iotools import ETSF_Reader
from abipy.abio.inputs import AnaddbInput
@ -65,6 +67,17 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
"""
This object provides an interface to the DDB file produced by ABINIT
as well as methods to compute phonon band structures, phonon DOS, thermodinamical properties ...
About the indices (idir, ipert) used by Abinit (Fortran notation)
* idir in [1, 2, 3] gives the direction (usually reduced direction)
* ipert in [1, 2, ..., mpert] where mpert = natom + 6
* ipert in [1, ..., natom] corresponds to atomic perturbations
* ipert = natom + 1 gives d/dk
* ipert = natom + 2 gives the electric field
* ipert = natom + 3 gives the uniaxial stress
* ipert = natom + 4 gives the shear stree.
"""
Error = DdbError
AnaddbError = AnaddbError
@ -102,10 +115,15 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
app("")
app(self.structure.to_string(verbose=verbose, title="Structure"))
app("")
app(self.qpoints.to_string(verbose=verbose, title="Q-points"))
app(self.qpoints.to_string(verbose=verbose, title="Q-points in DDB"))
app("")
app("guessed_ngqpt: %s (guess for the q-mesh divisions made by AbiPy)" % self.guessed_ngqpt)
if verbose > 1:
from pprint import pformat
app(marquee("DDB Header", mark="="))
app(pformat(self.header))
return "\n".join(lines)
@property
@ -160,14 +178,18 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
tokens = line.split()
key = None
try:
float(tokens[0])
parse = float if "." in tokens[0] else int
keyvals[-1][1].extend(list(map(parse, tokens)))
except ValueError:
# We have a new key
key = tokens.pop(0)
parse = float if "." in tokens[0] else int
keyvals.append((key, list(map(parse, tokens))))
try:
float(tokens[0])
parse = float if "." in tokens[0] else int
keyvals[-1][1].extend(list(map(parse, tokens)))
except ValueError:
# We have a new key
key = tokens.pop(0)
parse = float if "." in tokens[0] else int
keyvals.append((key, list(map(parse, tokens))))
except Exception as exc:
raise RuntimeError("Exception:\n%s\nwhile parsing ddb header line:\n%s" %
(str(exc), line))
# add the potential information
for line in self:
@ -184,17 +206,31 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
# to avoid problems with pymatgen routines that expect integral Z
# This of course will break any code for alchemical mixing.
arrays = {
"acell": dict(shape=(3, ), dtype=np.double),
"amu": dict(shape=(h.ntypat, ), dtype=np.double),
"kpt": dict(shape=(h.nkpt, 3), dtype=np.double),
"ngfft": dict(shape=(3, ), dtype=np.int),
# This is problematic because not all occupation factors are written
#"occ": dict(shape=(h.nsppol, h.nkpt, h.nband), dtype=np.double),
"rprim": dict(shape=(3, 3), dtype=np.double),
"spinat": dict(shape=(h.natom, 3), dtype=np.double),
"symrel": dict(shape=(h.nsym, 3, 3), dtype=np.int),
"tnons": dict(shape=(h.nsym, 3), dtype=np.double),
"xred": dict(shape=(h.natom, 3), dtype=np.double),
"znucl": dict(shape=(-1,), dtype=np.int),
"symafm": dict(shape=(-1,), dtype=np.int),
# In principle these two quantities are double but here we convert to int
# Alchemical mixing is therefore ignored.
"znucl": dict(shape=(h.ntypat,), dtype=np.int),
"zion": dict(shape=(h.ntypat,), dtype=np.int),
"symafm": dict(shape=(h.nsym,), dtype=np.int),
"wtk": dict(shape=(h.nkpt,), dtype=np.double),
}
for k, ainfo in arrays.items():
h[k] = np.reshape(np.array(h[k], dtype=ainfo["dtype"]), ainfo["shape"])
try:
h[k] = np.reshape(np.array(h[k], dtype=ainfo["dtype"]), ainfo["shape"])
except Exception as exc:
print("While Trying to reshape", k)
raise exc
# Transpose symrel because Abinit write matrices by colums.
h.symrel = np.array([s.T for s in h.symrel])
@ -225,33 +261,77 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
return np.reshape(qpoints, (-1,3))
@lazy_property
def computed_dynmat(self):
"""
OrderedDict mapping q-point object to --> pandas Dataframe.
The dataframe contains the columns: "idir1", "ipert1", "idir2", "ipert2", "cvalue"
and (idir1, ipert1, idir2, ipert2) as index.
.. note:
The indices follow the Abinit (Fortran) notation so they start at 1.
"""
# TODO: Create mapping [(idir1, ipert1), (idir2, ipert2)] --> element
df_columns = "idir1 ipert1 idir2 ipert2 cvalue".split()
dynmat = OrderedDict()
for block in self.blocks:
# Build q-point object.
qpt = Kpoint(frac_coords=block["qpt"], lattice=self.structure.reciprocal_lattice, weight=None, name=None)
# Build pandas dataframe with df_columns and (idir1, ipert1, idir2, ipert2) as index.
# Each line in data represents an element of the dynamical matric
# idir1 ipert1 idir2 ipert2 re_D im_D
df_rows, df_index = [], []
for line in block["data"]:
line = line.strip()
if line.startswith("2nd derivatives") or line.startswith("qpt"):
continue
try:
toks = line.split()
idir1, ipert1 = p1 = (int(toks[0]), int(toks[1]))
idir2, ipert2 = p2 = (int(toks[2]), int(toks[3]))
toks[4] = toks[4].replace("D", "E")
toks[5] = toks[5].replace("D", "E")
cvalue = float(toks[4]) + 1j*float(toks[5])
except Exception as exc:
print("exception while parsing line:", line)
raise exc
df_index.append(p1 + p2)
df_rows.append(dict(idir1=idir1, ipert1=ipert1, idir2=idir2, ipert2=ipert2, cvalue=cvalue))
dynmat[qpt] = pd.DataFrame(df_rows, index=df_index, columns=df_columns)
return dynmat
@lazy_property
def blocks(self):
"""
DDB blocks. List of dictionaries, Each dictionary contains "qpt"
with the reduced coordinates of the q-point and "data" that is a list of strings
with the entries of the dynamical matrix for this q-point.
DDB blocks. List of dictionaries, Each dictionary contains the following keys.
"qpt" with the reduced coordinates of the q-point.
"data" that is a list of strings with the entries of the dynamical matrix for this q-point.
"""
return self._read_blocks()
def _read_blocks(self):
self.seek(0)
# skip until the beginning of the db
self.seek(0)
while "Number of data blocks" not in self._file.readline():
pass
blocks = []
block_lines = []
qpt = None
for line in self:
# skip empty lines
if line.isspace():
continue
if "List of bloks and their characteristics" in line:
#add last block
# add last block when we reach the last part of the file.
blocks.append({"data": block_lines, "qpt": qpt})
break
@ -267,8 +347,6 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
if "qpt" in line:
qpt = list(map(float, line.split()[1:4]))
# TODO: Create mapping [(idir1, ipert1), (idir2, ipert2)] --> element
return blocks
@property
@ -276,6 +354,15 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
""":class:`KpointList` object with the list of q-points in reduced coordinates."""
return self._qpoints
def has_qpoint(self, qpoint):
"""True if the DDB file contains this q-point."""
#qpoint = Kpoint.as_kpoint(qpoint, self.structure.reciprocal_lattice)
try:
self.qpoints.index(qpoint)
return True
except ValueError:
return False
def qindex(self, qpoint):
"""
The index of the q-point in the internal list of q-points.
@ -331,37 +418,71 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
"""Dictionary with the parameters that are usually tested for convergence."""
return {k: v for k, v in self.header.items() if k in ("nkpt", "nsppol", "ecut", "tsmear", "ixc")}
# TODO
# API to understand if the DDB contains the info we are looking for.
# NB: This requires the parsing of the dynamical matrix
#def has_phonon_terms(self, qpoint)
# """True if the DDB file contains info on the phonon perturbation."""
def has_lo_to_data(self, select="at_least_one"):
"""
True if the DDB file contains data requires to compute LO-TO splitting.
"""
return self.has_emacro_terms(select=select) and self.has_bec_terms(select=select)
#def has_emacro_terms(self, ij="at_least_one"):
# """
# True if the DDB file contains info on the electric-field perturbation.
def has_emacro_terms(self, select="at_least_one"):
"""
True if the DDB file contains info on the electric-field perturbation.
# Args:
# ij: Possible values in ["at_least_one", "all"] or tuple e.g. (0, 1)
# If ij == "at_least_one", we check if there's at least one entry associated to the electric field.
# and we assume that anaddb will be able to reconstruct the full tensor by symmetry.
# If ij == "all", all tensor components must be present in the DDB file.
# If ij == (0, 1), the method returns False if the (0, 1) component of the tensor is not present in the DDB.
# """
Args:
select: Possible values in ["at_least_one", "all"]
If select == "at_least_one", we check if there's at least one entry associated to the electric field.
and we assume that anaddb will be able to reconstruct the full tensor by symmetry.
If select == "all", all tensor components must be present in the DDB file.
"""
gamma = Kpoint.gamma(self.structure.reciprocal_lattice)
if gamma not in self.computed_dynmat:
return False
#def has_bec_terms(self, ij="at_least_one"):
# """
# True if the DDB file contains info on the Born effective charges.
# By default, we check if there's at least one entry associated to electric field.
# and we assume that anaddb will be able to reconstruct the full tensor by symmetry
index_set = set(self.computed_dynmat[gamma].index)
# Args:
# ij: "at_least_one", "all", (1,2)
# """
natom = len(self.structure)
ep_list = list(itertools.product(range(1, 4), [natom + 2]))
for p1 in ep_list:
for p2 in ep_list:
p12 = p1 + p2
if select == "at_least_one":
if p12 in index_set: return True
elif select == "all":
if p12 not in index_set: return False
else:
raise ValueError("Wrong select %s" % str(select))
#def has_lo_to_data(self)
# """True if the DDB file contains data requires to compute LO-TO splitting."""
# return self.has_bec_terms() and self.has_emacro_terms()
return False
def has_bec_terms(self, select="at_least_one"):
"""
True if the DDB file contains info on the Born effective charges.
Args:
select: Possible values in ["at_least_one", "all"]
By default, we check if there's at least one entry associated to atomic displacement
and electric field and we assume that anaddb will be able to reconstruct the full tensor by symmetry.
If select == "all", all bec components must be present in the DDB file.
"""
gamma = Kpoint.gamma(self.structure.reciprocal_lattice)
if gamma not in self.computed_dynmat:
return False
index_set = set(self.computed_dynmat[gamma].index)
natom = len(self.structure)
ep_list = list(itertools.product(range(1, 4), [natom + 2]))
ap_list = list(itertools.product(range(1, 4), range(1, natom + 1)))
for ap1 in ap_list:
for ep2 in ep_list:
p12 = ap1 + ep2
if select == "at_least_one":
if p12 in index_set: return True
elif select == "all":
if p12 not in index_set: return False
else:
raise ValueError("Wrong select %s" % str(select))
return False
def anaget_phmodes_at_qpoint(self, qpoint=None, asr=2, chneut=1, dipdip=1, workdir=None, mpi_procs=1,
manager=None, verbose=0, lo_to_splitting=False, directions=None, anaddb_kwargs=None):
@ -397,7 +518,7 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
raise ValueError("input qpoint %s not in %s.\nddb.qpoints:\n%s" % (
qpoint, self.filepath, self.qpoints))
#if lo_to_splitting and not self.has_lo_to_data:
#if lo_to_splitting and qpoint.is_gamma() and not self.has_lo_to_data():
# cprint("lo_to_splitting set to True but Emacro and Becs are not available in DDB:" % self.filepath)
inp = AnaddbInput.modes_at_qpoint(self.structure, qpoint, asr=asr, chneut=chneut, dipdip=dipdip,
@ -422,12 +543,6 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
return ncfile.phbands
#def anaget_phbst_file(self, ngqpt=None, ndivsm=20, asr=2, chneut=1, dipdip=1,
# workdir=None, manager=None, verbose=0, **kwargs):
#def anaget_phdos_file(self, ngqpt=None, nqsmall=10, asr=2, chneut=1, dipdip=1, dos_method="tetra"
# workdir=None, manager=None, verbose=0, **kwargs):
def anaget_phbst_and_phdos_files(self, nqsmall=10, ndivsm=20, asr=2, chneut=1, dipdip=1, dos_method="tetra",
lo_to_splitting=False, ngqpt=None, qptbounds=None, anaddb_kwargs=None, verbose=0,
mpi_procs=1, workdir=None, manager=None):
@ -457,8 +572,8 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
"""
if ngqpt is None: ngqpt = self.guessed_ngqpt
#if lo_to_splitting and not self.has_lo_to_data:
# cprint("lo_to_splitting set to True but Emacro and Becs are not available in DDB:" % self.filepath)
if lo_to_splitting and not self.has_lo_to_data():
cprint("lo_to_splitting set to True but Emacro and Becs are not available in DDB:" % self.filepath, "yellow")
inp = AnaddbInput.phbands_and_dos(
self.structure, ngqpt=ngqpt, ndivsm=ndivsm, nqsmall=nqsmall, q1shft=(0, 0, 0), qptbounds=qptbounds,
@ -511,7 +626,6 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
if num_cpus <= 0: num_cpus = 1
num_cpus = min(num_cpus, len(nqsmalls))
# TODO: anaget_phdos
def do_work(nqsmall):
_, phdos_file = self.anaget_phbst_and_phdos_files(
nqsmall=nqsmall, ndivsm=1, asr=asr, chneut=chneut, dipdip=dipdip, dos_method=dos_method, ngqpt=ngqpt)
@ -580,6 +694,9 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
Return:
emacro, becs
"""
if not self.has_lo_to_data():
cprint("Dielectric tensor and Becs are not available in DDB: %s" % self.filepath, "yellow")
inp = AnaddbInput(self.structure, anaddb_kwargs={"chneut": chneut})
task = AnaddbTask.temp_shell_task(inp, ddb_node=self.filepath, mpi_procs=mpi_procs, workdir=workdir, manager=manager)
@ -598,6 +715,7 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
with ETSF_Reader(os.path.join(task.workdir, "anaddb.nc")) as r:
structure = r.read_structure()
# TODO Replace with pymatgen tensors
emacro = Tensor.from_cartesian_tensor(r.read_value("emacro_cart"), structure.lattice, space="r"),
becs = Becs(r.read_value("becs_cart"), structure, chneut=inp["chneut"], order="f")
@ -673,10 +791,10 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
break
header_lines.append(line)
nskip=10
nskip = 10
fmti = "%5d"
fmtf = "%22.14e"
fmt3 = " "*nskip+fmtf*3+'\n'
fmt3 = " " * nskip + fmtf *3 + '\n'
#write all the variables in order
for variable in variables:
@ -1143,7 +1261,7 @@ class DielectricTensorGenerator(Has_Structure):
return DielectricTensor(t)
@add_fig_kwargs
def plot_vs_w(self, w_min, w_max, num, component='diago', units='eV', ax=None, **kwargs):
def plot_vs_w(self, w_min, w_max, num, component='diago', units='eV', ax=None, fontsize=12, **kwargs):
"""
Plots the selected components of the dielectric tensor as a function of the frequency
@ -1153,10 +1271,11 @@ class DielectricTensorGenerator(Has_Structure):
num: number of values of the frequencies between w_min and w_max
component: determine which components of the tensor will be displayed. Can be a list/tuple of two
elements, indicating the indices [i, j] of the desired component or a string among
'diag': plots the elements on diagonal
'all': plots all the components
'diag_av': plots the average of the components on the diagonal
'diag' to plot the elements on diagonal
'all' to plot all the components
'diag_av' to plot the average of the components on the diagonal
units: string specifying the units used for the frequency. Accepted values are Ha, eV (default), cm-1
fontsize: Legend and label fontsize.
"""
w_range = np.linspace(w_min, w_max, num, endpoint=True)
@ -1194,7 +1313,7 @@ class DielectricTensorGenerator(Has_Structure):
else:
ValueError('Unkwnown component {}'.format(component))
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@ -1249,11 +1368,11 @@ class DdbRobot(Robot):
"""
# If qpoint is None, all the DDB must contain have the same q-point .
if qpoint is None:
if not all(len(ddb.qpoints) == 1 for ddb in self.ncfiles):
if not all(len(ddb.qpoints) == 1 for ddb in self.abifiles):
raise ValueError("Found more than one q-point in the DDB file. qpoint must be specified")
qpoint = self[0].qpoints[0]
if any(np.any(ddb.qpoints[0] != qpoint) for ddb in self.ncfiles):
if any(np.any(ddb.qpoints[0] != qpoint) for ddb in self.abifiles):
raise ValueError("All the q-points in the DDB files must be equal")
rows, row_names = [], []
@ -1280,31 +1399,10 @@ class DdbRobot(Robot):
if funcs is not None: d.update(self._exec_funcs(funcs, ddb))
rows.append(d)
import pandas as pd
row_names = row_names if not abspath else self._to_relpaths(row_names)
return pd.DataFrame(rows, index=row_names, columns=list(rows[0].keys()))
# TODO: Is this really needed?
#def plot_conv_phfreqs_qpoint(self, x_vars, qpoint=None, **kwargs):
# """
# Plot the convergence of the phonon frequencies.
# kwargs are passed to :class:`seaborn.PairGrid`.
# """
# import matplotlib.pyplot as plt
# import seaborn.apionly as sns
# # Get the dataframe for this q-point.
# data = self.get_dataframe_at_qpoint(qpoint=qpoint)
# y_vars = sorted([k for k in data if k.startswith("mode")])
# #print(y_vars)
# # Call seaborn.
# grid = sns.PairGrid(data, x_vars=x_vars, y_vars=y_vars, **kwargs)
# grid.map(plt.plot, marker="o")
# grid.add_legend()
# plt.show()
def anaget_phonon_plotters(self, **kwargs):
r"""
Invoke anaddb to compute phonon bands and DOS using the arguments passed via \*\*kwargs.

2
abipy/dfpt/gruneisen.py
View File

@ -180,7 +180,7 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
if i == len(dos_names) - 1:
ax.set_xlabel(r"$\omega$ [eV]")
#ax.legend(loc="best")
#ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig

175
abipy/dfpt/phonons.py
View File

@ -437,10 +437,16 @@ class PhononBands(object):
Write xmgrace file with phonon band structure energies and labels for high-symmetry q-points.
Args:
filepath: Filename
filepath: String with filename or stream.
units: Units for phonon plots. Possible values in ("eV", "meV", "Ha", "cm-1", "Thz"). Case-insensitive.
"""
f = open(filepath, "wt")
is_stream = hasattr(filepath, "write")
if is_stream:
f = filepath
else:
f = open(filepath, "wt")
def w(s):
f.write(s)
f.write("\n")
@ -498,7 +504,8 @@ class PhononBands(object):
w('%d %.8E' % (iq, wqnu_units[iq, nu]))
w('&')
f.close()
if not is_stream:
f.close()
#def to_bxsf(self, filepath):
# """
@ -711,43 +718,29 @@ class PhononBands(object):
with open(filename, 'wt') as f:
f.write("\n".join(lines))
def view_phononwebsite(self, open_browser=True, verbose=1, **kwargs):
def view_phononwebsite(self, browser=None, verbose=1, **kwargs):
"""
Produce json file that can be parsed from the phononwebsite. Contact the server, get
the url of the webpage and open it in the browser if `open_browser`.
Produce json file that can be parsed from the phononwebsite and open it in `browser`.
Return:
Exit status
"""
filename = kwargs.pop("filename", None)
if filename is None:
import tempfile
prefix = self.structure.formula.replace(" ", "")
_, filename = tempfile.mkstemp(text=True, prefix=prefix, suffix=".json")
# Create json in abipy_nbworkdir with relative path so that we can read it inside the browser.
from abipy.core.globals import abinb_mkstemp
prefix = self.structure.formula.replace(" ", "")
_, rpath = abinb_mkstemp(force_abinb_workdir=True, use_relpath=True,
prefix=prefix, suffix=".json", text=True)
if verbose: print("Writing json file", filename)
self.create_phononwebsite_json(filename, indent=None, **kwargs)
url = "http://henriquemiranda.github.io/phononwebsite/phonon.html"
import requests
with open(filename, 'rt') as f:
files = {'file-input': (filename, f)}
r = requests.post(url, files=files)
# @Henrique: the response should contain the url of the bandstructure plot
# so that I can open it in the browser.
if verbose:
print(r)
print(r.text)
if verbose: print("Writing json file:", rpath)
self.create_phononwebsite_json(rpath, indent=None, **kwargs)
print("Phonon band structure available at:", phbst_url)
if open_browser:
import webbrowser
return int(webbrowser.open(phbst_url))
return 0
return open_file_phononwebsite(rpath, browser=browser)
def create_phononwebsite_json(self, filename, name=None, repetitions=None, highsym_qpts=None, match_bands=True,
highsym_qpts_mode="split", indent=2):
def create_phononwebsite_json(self, filename, name=None, repetitions=None, highsym_qpts=None,
match_bands=True, highsym_qpts_mode="split", indent=2):
"""
Writes a json file that can be parsed from the phononwebsite. See <https://github.com/henriquemiranda/phononwebsite>
Writes a json file that can be parsed from the phononwebsite.
See <https://github.com/henriquemiranda/phononwebsite>
Args:
filename: name of the json file that will be created
@ -954,7 +947,8 @@ class PhononBands(object):
units: Units for phonon plots. Possible values in ("eV", "meV", "Ha", "cm-1", "Thz"). Case-insensitive.
qlabels: dictionary whose keys are tuples with the reduced coordinates of the q-points.
The values are the labels. e.g. ``qlabels = {(0.0,0.0,0.0): "$\Gamma$", (0.5,0,0): "L"}``
branch_range: Tuple specifying the minimum and maximum branch_i index to plot (default: all branches are plotted).
branch_range: Tuple specifying the minimum and maximum branch_i index to plot
(default: all branches are plotted).
colormap: matplotlib colormap to determine the colors available. The colors will be chosen not in a
sequential order to avoid difficulties in distinguishing the lines.
http://matplotlib.sourceforge.net/examples/pylab_examples/show_colormaps.html
@ -1186,7 +1180,7 @@ class PhononBands(object):
@add_fig_kwargs
def plot_fatbands(self, units="eV", colormap="jet", phdos_file=None,
alpha=0.7, max_stripe_width_mev=3.0, width_ratios=(2, 1),
qlabels=None, ylims=None,
qlabels=None, ylims=None, fontsize=12,
**kwargs):
#cart_dir=None
r"""
@ -1206,6 +1200,7 @@ class PhononBands(object):
or scalar e.g. `left`. If left (right) is None, default values are used
qlabels: dictionary whose keys are tuples with the reduced coordinates of the q-points.
The values are the labels. e.g. ``qlabels = {(0.0,0.0,0.0): "$\Gamma$", (0.5,0,0): "L"}``.
fontsize: Legend and title fontsize.
Returns:
`matplotlib` figure.
@ -1287,7 +1282,7 @@ class PhononBands(object):
ax.fill_between(qq, yy_qq + d2_type, yy_qq - d2_type, facecolor=color, alpha=alpha, linewidth=0)
set_axlims(ax, ylims, "y")
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
# Type projected DOSes.
if phdos_file is not None:
@ -1923,7 +1918,7 @@ class PhononDos(Function1D):
ax.grid(True)
ax.set_xlabel("Temperature [K]")
ax.set_ylabel(_THERMO_YLABELS[qname][units])
#ax.legend(loc="best")
#ax.legend(loc="best", fontsize=fontsize, shadow=True)
fig.tight_layout()
return fig
@ -2065,7 +2060,7 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
@add_fig_kwargs
def plot_pjdos_type(self, units="eV", stacked=True, colormap="jet", alpha=0.7,
ax=None, xlims=None, ylims=None, **kwargs):
ax=None, xlims=None, ylims=None, fontsize=12, **kwargs):
"""
Plot type-projected phonon DOS.
@ -2080,6 +2075,7 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
xlims: Set the data limits for the x-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
ylims: y-axis limits.
fontsize: legend and title fontsize.
Returns:
matplotlib figure.
@ -2113,13 +2109,13 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
# Total PHDOS
x, y = self.phdos.mesh * factor, self.phdos.values / factor
ax.plot(x, y, lw=lw, label="Total PHDOS", color='black')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@add_fig_kwargs
def plot_pjdos_redirs_type(self, units="eV", stacked=True, colormap="jet", alpha=0.7,
xlims=None, ylims=None, axlist=None, **kwargs):
xlims=None, ylims=None, axlist=None, fontsize=12, **kwargs):
"""
Plot type-projected phonon DOS decomposed along the three reduced directions.
Three rows for each reduced direction. Each row shows the contribution of each atomic type + Total PH DOS.
@ -2135,6 +2131,7 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
or scalar e.g. `left`. If left (right) is None, default values are used
ylims: y-axis limits.
axlist: List of matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: Legend and label fontsize.
Returns:
matplotlib figure.
@ -2182,13 +2179,13 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
# Add Total PHDOS
ax.plot(xx, self.phdos.values / factor, lw=lw, label="Total PHDOS", color='black')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@add_fig_kwargs
def plot_pjdos_redirs_site(self, view="inequivalent", units="eV", stacked=True, colormap="jet", alpha=0.7,
xlims=None, ylims=None, axlist=None, **kwargs):
xlims=None, ylims=None, axlist=None, fontsize=12, **kwargs):
"""
Plot phonon PJDOS for each atom in the unit cell. By default, only "inequivalent" atoms are shown.
@ -2203,6 +2200,7 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
ylims: Set the data limits for the y-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
axlist: List of matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: Legend and title fontsize.
Returns:
`matplotlib` figure
@ -2264,7 +2262,7 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
# Add Total PHDOS
ax.plot(xx, self.phdos.values / factor, lw=lw, label="Total PHDOS", color='black')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@ -2879,7 +2877,7 @@ class PhononDosPlotter(NotebookWriter):
self._phdoses_dict[label] = PhononDos.as_phdos(phdos, phdos_kwargs)
@add_fig_kwargs
def combiplot(self, ax=None, units="eV", xlims=None, ylims=None, **kwargs):
def combiplot(self, ax=None, units="eV", xlims=None, ylims=None, fontsize=12, **kwargs):
"""
Plot DOSes on the same figure. Use `gridplot` to plot DOSes on different figures.
@ -2889,6 +2887,7 @@ class PhononDosPlotter(NotebookWriter):
xlims: Set the data limits for the x-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
ylims: y-axis limits.
fontsize: Legend and title fontsize.
"""
ax, fig, plt = get_ax_fig_plt(ax)
@ -2905,7 +2904,7 @@ class PhononDosPlotter(NotebookWriter):
legends.append("DOS: %s" % label)
# Set legends.
ax.legend(lines, legends, loc='best', shadow=True)
ax.legend(lines, legends, loc='best', fontsize=fontsize, shadow=True)
return fig
@ -2959,7 +2958,7 @@ class PhononDosPlotter(NotebookWriter):
@add_fig_kwargs
def plot_harmonic_thermo(self, tstart=5, tstop=300, num=50, units="eV", formula_units=1,
quantities="all", **kwargs):
quantities="all", fontsize=12, **kwargs):
"""
Plot thermodinamic properties from the phonon DOS within the harmonic approximation.
@ -2973,6 +2972,7 @@ class PhononDosPlotter(NotebookWriter):
thermodynamic quantities will be given on a per-unit-cell basis.
quantities: List of strings specifying the thermodinamic quantities to plot.
Possible values: ["internal_energy", "free_energy", "entropy", "c_v"].
fontsize: Legend and title fontsize.
Returns:
matplotlib figure.
@ -3002,11 +3002,11 @@ class PhononDosPlotter(NotebookWriter):
if units == "Jmol": ys = ys * abu.e_Cb * abu.Avogadro
ax.plot(f1d.mesh, ys, label=label)
ax.set_title(qname)
ax.set_title(qname, fontsize=fontsize)
ax.grid(True)
ax.set_xlabel("Temperature [K]")
ax.set_ylabel(_THERMO_YLABELS[qname][units])
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
fig.tight_layout()
return fig
@ -3513,7 +3513,7 @@ class RobotWithPhbands(object):
"""
Build a pandas dataframe with the most important results available in the band structures.
"""
return dataframe_from_phbands([nc.phbands for nc in self.ncfiles],
return dataframe_from_phbands([nc.phbands for nc in self.abifiles],
index=self.labels, with_spglib=with_spglib)
def get_phbands_code_cells(self, title=None):
@ -3525,3 +3525,84 @@ class RobotWithPhbands(object):
nbv.new_markdown_cell(title),
nbv.new_code_cell("robot.get_phbands_plotter().ipw_select_plot();"),
]
def open_file_phononwebsite(filename, port=8000,
website="http://henriquemiranda.github.io/phononwebsite",
host="localhost", browser=None):
"""
Take a file, detect the type and open it on the phonon website
Based on a similar function in https://github.com/henriquemiranda/phononwebsite/phononweb.py
Args:
Return
"""
if filename.endswith(".json"):
filetype = "json"
elif filename.endswith(".yaml"):
filetype = "yaml"
else:
filetype = "rest"
try:
from http.server import HTTPServer, SimpleHTTPRequestHandler
except ImportError:
from BaseHTTPServer import HTTPServer
# python 2 requires internal implementation
from abipy.tools.SimpleHTTPServer import SimpleHTTPRequestHandler
# Add CORS header to the website
class CORSRequestHandler (SimpleHTTPRequestHandler):
def end_headers (self):
#self.send_header('Access-Control-Allow-Origin', website)
self.send_header('Access-Control-Allow-Origin', "http://henriquemiranda.github.io")
SimpleHTTPRequestHandler.end_headers(self)
def log_message(self, format, *args):
return
# Initialize http server thread
print('Starting HTTP server at port %d ...' % port, end=" ")
trial, max_ntrial = 0, 50
while trial < max_ntrial:
try:
server = HTTPServer(('', port), CORSRequestHandler)
#print("got port:", port)
break
except OSError:
trial += 1
port += 1
print(port, end=", ")
else:
raise RuntimeError("Cannot find available port after %s attempts" % max_ntrial)
# Create threads python
server.url = 'http://{}:{}'.format(host, server.server_port)
from threading import Thread
t = Thread(target=server.serve_forever)
t.daemon = True
t.start()
# Open website with the file
try:
from urllib.parse import quote
except ImportError:
from urllib import quote
url_filename = 'http://{}:{}/{}'.format(host, server.server_port, quote(filename))
#url_filename = 'http://{}:{}/{}'.format(host, server.server_port, filename)
url = '%s/phonon.html?%s=%s' % (website, filetype, url_filename)
print("\nOpening URL:", url)
print("Using default browser, if the webpage is not displayed correctly",
"\ntry to change browser either via command line options or directly in the shell with e.g:\n\n"
" export BROWSER=firefox\n")
print('Press Ctrl+C to terminate HTTP server')
import webbrowser
webbrowser.get(browser).open_new(url)
# Quit application when SIGINT is received
import sys
def signal_handler(signal, frame):
sys.exit(0)
import signal
signal.signal(signal.SIGINT, signal_handler)
signal.pause()

2
abipy/dfpt/tests/test_anaddbnc.py
View File

@ -17,7 +17,7 @@ class AnaddbNcFileTest(AbipyTest):
with AnaddbNcFile(anaddbnc_fname) as anc:
repr(anc); str(anc)
anc.to_string(verbose=1)
anc.to_string(verbose=2)
assert anc.structure.formula == "Al1 As1"
assert anc.becs is not None
assert anc.emacro is not None

81
abipy/dfpt/tests/test_ddb.py
View File

@ -45,6 +45,9 @@ class DdbTest(AbipyTest):
assert struct.formula == "Al1 As1"
# Test update header
# FIXME Disabled after my changes in ddb.header
# The implementation of update_header should be checked again.
"""
h_copy = ddb.header.copy()
ddb.update_header()
for k, v in ddb.header.items():
@ -56,6 +59,7 @@ class DdbTest(AbipyTest):
# err += 1
# print("line1", line1, "\nline2, line2)
#assert err == 0
"""
# Test interface with Anaddb.
print(ddb.qpoints[0])
@ -114,12 +118,29 @@ class DdbTest(AbipyTest):
"""Testing DDB for AlAs on a 4x4x4x q-mesh without Born effective charges."""
ddb = DdbFile(os.path.join(test_dir, "AlAs_444_nobecs_DDB"))
repr(ddb); str(ddb)
print(ddb.header)
assert str(ddb.header)
assert ddb.to_string(verbose=2)
assert ddb.header["nkpt"] == 256
assert ddb.header.nsym == 24 and ddb.header.ntypat == 2
self.assert_equal(ddb.header.znucl, [13, 33])
self.assert_equal(ddb.header.acell, [1, 1, 1])
self.assert_equal(ddb.header.ngfft, [10, 10, 10])
self.assert_equal(ddb.header.spinat, 0.0)
#assert ddb.header.occ.shape = (ddb.header.nsppol, ddb.header.nkpt, ddb.header.nsppol)
# TODO
#assert ddb.has_phonon_terms()
#assert not ddb.has_bec_terms()
#assert not ddb.has_emacro_terms()
assert not ddb.has_qpoint([0.345, 0.456, 0.567])
assert ddb.has_qpoint([0, 0, 0])
for qpoint in ddb.qpoints:
assert ddb.has_qpoint(qpoint)
assert ddb.has_qpoint(qpoint.frac_coords)
assert qpoint in ddb.computed_dynmat
assert len(ddb.computed_dynmat[qpoint].index[0]) == 4
assert ddb.has_bec_terms(select="at_least_one")
assert not ddb.has_bec_terms(select="all")
assert not ddb.has_emacro_terms()
assert not ddb.has_lo_to_data()
ref_qpoints = np.reshape([
0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
@ -153,9 +174,9 @@ class DdbTest(AbipyTest):
# Get emacro and becs
emacro, becs = ddb.anaget_emacro_and_becs(chneut=1, verbose=1)
assert np.all(becs.values == 0)
assert np.all(becs.values == 0)
#assert np.all(emacro.values == 0)
repr(becs); str(becs)
assert becs.to_string(verbose=1)
assert becs.to_string(verbose=2)
self.assert_almost_equal(phdos.idos.values[-1], 3 * len(ddb.structure), decimal=1)
phbands_file.close()
@ -180,6 +201,54 @@ class DdbTest(AbipyTest):
ddb.close()
def test_zno_gamma_ddb_with_becs(self):
"""Testing DDB for ZnO: Gamma only, with Born effective charges and E_macro."""
with DdbFile(os.path.join(test_dir, "ZnO_gamma_becs_DDB")) as ddb:
repr(ddb); str(ddb)
assert str(ddb.header)
assert ddb.to_string(verbose=2)
assert ddb.header["nkpt"] == 486
assert ddb.header.nband == 22 and ddb.header.occopt == 1
self.assert_equal(ddb.header.typat, [1, 1, 2, 2])
assert len(ddb.header.wtk) == ddb.header.nkpt
#assert ddb.header.occ.shape = (ddb.header.nsppol, ddb.header.nkpt, ddb.header.nsppol)
assert not ddb.has_qpoint([0.345, 0.456, 0.567])
assert ddb.has_qpoint([0, 0, 0])
assert len(ddb.qpoints) == 1
for qpoint in ddb.qpoints:
assert ddb.has_qpoint(qpoint)
assert ddb.has_qpoint(qpoint.frac_coords)
assert qpoint in ddb.computed_dynmat
assert len(ddb.computed_dynmat[qpoint].index[0]) == 4
assert ddb.has_bec_terms(select="at_least_one")
assert not ddb.has_bec_terms(select="all")
assert ddb.has_emacro_terms()
assert ddb.has_lo_to_data()
# Get emacro and becs
emacro, becs = ddb.anaget_emacro_and_becs(chneut=1, verbose=1)
ref_becs_values = [
[[ 2.15646571e+00, 0.00000000e+00, 3.26402110e-25],
[ 0.00000000e+00, 2.15646571e+00, -5.46500204e-24],
[ -5.66391495e-25, -6.54012564e-25, 2.19362823e+00]],
[[ 2.15646571e+00, 0.00000000e+00, 1.19680774e-24],
[ 0.00000000e+00, 2.15646571e+00, 8.10327888e-24],
[ -1.69917448e-24, -1.30802513e-24, 2.19362823e+00]],
[[ -2.15646571e+00, 6.66133815e-16, -1.84961196e-24],
[ 8.88178420e-16, -2.15646571e+00, 2.82672519e-24],
[ -3.39834897e-24, -3.27006282e-25, -2.19362823e+00]],
[[ -2.15646571e+00, -6.66133815e-16, 3.26402110e-25],
[ -8.88178420e-16, -2.15646571e+00, -5.46500204e-24],
[ 5.66391495e-24, 2.28904397e-24, -2.19362823e+00]]
]
self.assert_almost_equal(becs.values, ref_becs_values)
#self.assert_almost_equal(emacro.values, ref_emacro_values)
repr(becs); str(becs)
assert becs.to_string(verbose=2)
class DielectricTensorGeneratorTest(AbipyTest):

2
abipy/dfpt/tests/test_phonons.py
View File

@ -2,6 +2,7 @@
from __future__ import print_function, division, unicode_literals, absolute_import
import unittest
import sys
import os
import pickle
import numpy as np
@ -78,6 +79,7 @@ class PhononBandsTest(AbipyTest):
phbands.create_phononwebsite_json(filename=self.get_tmpname(text=True), name='test')
# Test xmgrace
phbands.to_xmgrace(self.get_tmpname(text=True))
phbands.to_xmgrace(sys.stdout)
df = phbands.get_dataframe()
assert "freq" in df and "mode" in df

340
abipy/dynamics/hist.py
View File

@ -9,9 +9,9 @@ import pymatgen.core.units as units
from collections import OrderedDict
from monty.functools import lazy_property
from monty.collections import AttrDict
from monty.string import marquee # is_string, list_strings
from monty.string import marquee, list_strings
from pymatgen.core.periodic_table import Element
from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt
from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, get_axarray_fig_plt
from abipy.core.structure import Structure
from abipy.core.mixins import AbinitNcFile, NotebookWriter
from abipy.abio.robots import Robot
@ -44,6 +44,23 @@ class HistFile(AbinitNcFile, NotebookWriter):
def __str__(self):
return self.to_string()
# TODO
#@lazy_property
#def nsppol(self):
# """Number of independent spins."""
# return self.reader.read_dimvalue("nsppol")
#@lazy_property
#def nspden(self):
# """Number of independent spin densities."""
# return self.reader.read_dimvalue("nspden")
#@lazy_property
#def nspinor(self):
# """Number of spinor components."""
# return self.reader.read_dimvalue("nspinor")
@lazy_property
def final_energy(self):
return self.etotals[-1]
@ -57,17 +74,22 @@ class HistFile(AbinitNcFile, NotebookWriter):
#@lazy_property
#def final_max_force(self):
#def get_fstats_dict(self, step):
# #for step in range(self.num_steps):
# forces_hist = self.reader.read_cart_forces()
# fmin_steps, fmax_steps, fmean_steps, fstd_steps = [], [], [], []
def get_fstats_dict(self, step):
"""
Return dictionary with stats on the forces at the given `step`
"""
# [time, natom, 3]
var = self.reader.read_variable("fcart")
forces = units.ArrayWithUnit(var[step], "Ha bohr^-1").to("eV ang^-1")
fmods = np.array([np.linalg.norm(force) for force in forces])
# forces = forces_hist[step]
# fmods = np.sqrt([np.dot(force, force) for force in forces])
# fmean_steps.append(fmods.mean())
# fstd_steps.append(fmods.std())
# fmin_steps.append(fmods.min())
# fmax_steps.append(fmods.max())
return AttrDict(
fmin=fmods.min(),
fmax=fmods.max(),
fmean=fmods.mean(),
fstd=fmods.std(),
drift=np.linalg.norm(forces.sum(axis=0)),
)
def to_string(self, verbose=0, title=None):
"""Return string representation."""
@ -164,13 +186,19 @@ class HistFile(AbinitNcFile, NotebookWriter):
raise RuntimeError("Cannot overwrite pre-existing file `%s`" % filepath)
if filepath is None:
import tempfile
fd, filepath = tempfile.mkstemp(text=True)
fd, filepath = tempfile.mkstemp(text=True, suffix="_XDATCAR")
# int typat[natom], double znucl[npsp]
typat = self.reader.read_value("typat")
# NB: typat is double in the HIST.nc file
typat = self.reader.read_value("typat").astype(int)
znucl = self.reader.read_value("znucl")
if len(typat) != len(znucl):
raise RuntimeError("Alchemical mixing is not supported.")
ntypat = self.reader.read_dimvalue("ntypat")
num_pseudos = self.reader.read_dimvalue("npsp")
if num_pseudos != ntypat:
raise NotImplementedError("Alchemical mixing is not supported, num_pseudos != ntypat")
print("znucl:", znucl)
print("typat:", typat)
symb2pos = OrderedDict()
symbols_atom = []
for iatom, itype in enumerate(typat):
@ -212,6 +240,119 @@ class HistFile(AbinitNcFile, NotebookWriter):
return filepath
def visualize(self, appname="ovito"):
"""
Visualize the crystalline structure with visualizer.
See :class:`Visualizer` for the list of applications and formats supported.
"""
if appname == "mayavi": return self.mayaview()
# Get the Visualizer subclass from the string.
from abipy.iotools import Visualizer
visu = Visualizer.from_name(appname)
if visu.name != "ovito":
raise NotImplementedError("visualizer: %s" % visu.name)
filepath = self.write_xdatcar(filepath=None, groupby_type=True)
return visu(filepath)()
#if options.trajectories:
# hist.mvplot_trajectories()
#else:
# hist.mvanimate()
def plot_ax(self, ax, what, fontsize=12, **kwargs):
"""
Helper function to plot quantity `what` on axis `ax`.
Args:
fontsize: fontsize for legend
kwargs are passed to matplotlib plot method
"""
if what == "energy":
# Total energy in eV.
marker = kwargs.pop("marker", "o")
label = kwargs.pop("label", "Energy")
ax.plot(self.steps, self.etotals, label=label, marker=marker, **kwargs)
ax.set_ylabel('Energy [eV]')
elif what == "abc":
# Lattice parameters.
mark = kwargs.pop("marker", None)
markers = ["o", "^", "v"] if mark is None else 3 * [mark]
for i, label in enumerate(["a", "b", "c"]):
ax.plot(self.steps, [s.lattice.abc[i] for s in self.structures], label=label,
marker=markers[i], **kwargs)
ax.set_ylabel("abc [A]")
elif what in ("a", "b", "c"):
i = ("a", "b", "c").index(what)
marker = kwargs.pop("marker", None)
if marker is None:
marker = {"a": "o", "b": "^", "c": "v"}[what]
label = kwargs.pop("label", what)
ax.plot(self.steps, [s.lattice.abc[i] for s in self.structures], label=label,
marker=marker, **kwargs)
ax.set_ylabel('%s [A]' % what)
elif what == "angles":
# Lattice Angles
mark = kwargs.pop("marker", None)
markers = ["o", "^", "v"] if mark is None else 3 * [mark]
for i, label in enumerate(["alpha", "beta", "gamma"]):
ax.plot(self.steps, [s.lattice.angles[i] for s in self.structures], label=label,
marker=markers[i], **kwargs)
ax.set_ylabel(r"$\alpha\beta\gamma$ [degree]")
elif what in ("alpha", "beta", "gamma"):
i = ("alpha", "beta", "gamma").index(what)
marker = kwargs.pop("marker", None)
if marker is None:
marker = {"alpha": "o", "beta": "^", "gamma": "v"}[what]
label = kwargs.pop("label", what)
ax.plot(self.steps, [s.lattice.angles[i] for s in self.structures], label=label,
marker=marker, **kwargs)
ax.set_ylabel(r"$\%s$ [degree]" % what)
elif what == "volume":
marker = kwargs.pop("marker", "o")
ax.plot(self.steps, [s.lattice.volume for s in self.structures], marker=marker, **kwargs)
ax.set_ylabel(r'$V\, [A^3]$')
elif what == "pressure":
stress_cart_tensors, pressures = self.reader.read_cart_stress_tensors()
marker = kwargs.pop("marker", "o")
label = kwargs.pop("label", "P")
ax.plot(self.steps, pressures, label=label, marker=marker, **kwargs)
ax.set_ylabel('P [GPa]')
elif what == "forces":
forces_hist = self.reader.read_cart_forces()
fmin_steps, fmax_steps, fmean_steps, fstd_steps = [], [], [], []
for step in range(self.num_steps):
forces = forces_hist[step]
fmods = np.sqrt([np.dot(force, force) for force in forces])
fmean_steps.append(fmods.mean())
fstd_steps.append(fmods.std())
fmin_steps.append(fmods.min())
fmax_steps.append(fmods.max())
mark = kwargs.pop("marker", None)
markers = ["o", "^", "v", "X"] if mark is None else 4 * [mark]
ax.plot(self.steps, fmin_steps, label="min |F|", marker=markers[0], **kwargs)
ax.plot(self.steps, fmax_steps, label="max |F|", marker=markers[1], **kwargs)
ax.plot(self.steps, fmean_steps, label="mean |F|", marker=markers[2], **kwargs)
ax.plot(self.steps, fstd_steps, label="std |F|", marker=markers[3], **kwargs)
ax.set_ylabel('F stats [eV/A]')
else:
raise ValueError("Invalid value for what: `%s`" % str(what))
ax.set_xlabel('Step')
ax.legend(loc='best', fontsize=fontsize, shadow=True)
ax.grid(True)
@add_fig_kwargs
def plot(self, axlist=None, **kwargs):
"""
@ -225,63 +366,24 @@ class HistFile(AbinitNcFile, NotebookWriter):
`matplotlib` figure
"""
import matplotlib.pyplot as plt
what_list = ["abc", "angles", "volume", "pressure", "forces", "energy"]
fig, ax_list = plt.subplots(nrows=3, ncols=2, sharex=True, squeeze=False)
ax_list = ax_list.ravel()
ax0, ax1, ax2, ax3, ax4, ax5 = ax_list
for ax in ax_list: ax.grid(True)
assert len(ax_list) == len(what_list)
# Lattice parameters.
for i, label in enumerate(["a", "b", "c"]):
ax0.plot(self.steps, [s.lattice.abc[i] for s in self.structures], marker="o", label=label)
ax0.set_ylabel('Lattice lengths [A]')
ax0.legend(loc='best', shadow=True)
# Lattice Angles
for i, label in enumerate(["alpha", "beta", "gamma"]):
ax1.plot(self.steps, [s.lattice.angles[i] for s in self.structures], marker="o", label=label)
ax1.set_ylabel('Lattice Angles [degree]')
ax1.legend(loc='best', shadow=True)
ax2.plot(self.steps, [s.lattice.volume for s in self.structures], marker="o")
ax2.set_ylabel('Lattice volume [A^3]')
stress_cart_tensors, pressures = self.reader.read_cart_stress_tensors()
ax3.plot(self.steps, pressures, marker="o", label="Pressure")
ax3.set_ylabel('Pressure [GPa]')
# Forces
forces_hist = self.reader.read_cart_forces()
fmin_steps, fmax_steps, fmean_steps, fstd_steps = [], [], [], []
for step in range(self.num_steps):
forces = forces_hist[step]
fmods = np.sqrt([np.dot(force, force) for force in forces])
fmean_steps.append(fmods.mean())
fstd_steps.append(fmods.std())
fmin_steps.append(fmods.min())
fmax_steps.append(fmods.max())
ax4.plot(self.steps, fmin_steps, marker="o", label="min |F|")
ax4.plot(self.steps, fmax_steps, marker="o", label="max |F|")
ax4.plot(self.steps, fmean_steps, marker="o", label="mean |F|")
ax4.plot(self.steps, fstd_steps, marker="o", label="std |F|")
ax4.set_ylabel('Force stats [eV/A]')
ax4.legend(loc='best', shadow=True)
ax4.set_xlabel('Step')
# Total energy.
ax5.plot(self.steps, self.etotals, marker="o", label="Energy")
ax5.set_ylabel('Total energy [eV]')
ax5.set_xlabel('Step')
for what, ax in zip(what_list, ax_list):
self.plot_ax(ax, what, marker="o")
return fig
@add_fig_kwargs
def plot_energies(self, ax=None, **kwargs):
def plot_energies(self, ax=None, fontsize=12, **kwargs):
"""
Plot the total energies as function of the iteration step.
Args:
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: Legend and title fontsize.
Returns:
`matplotlib` figure
@ -297,7 +399,7 @@ class HistFile(AbinitNcFile, NotebookWriter):
ax.set_xlabel('Step')
ax.set_ylabel('Energies [eV]')
ax.grid(True)
ax.legend(loc='best', shadow=True)
ax.legend(loc='best', fontsize=fontsize, shadow=True)
return fig
@ -397,6 +499,7 @@ class HistRobot(Robot):
EXT = "HIST"
def to_string(self, verbose=0):
"""String representation with verbosity level `verbose`."""
s = ""
if verbose:
s = super(HistRobot, self).to_string(verbose=0)
@ -427,11 +530,13 @@ class HistRobot(Robot):
where key is a string with the name of column and value is the value to be inserted.
"""
# Add attributes specified by the users
# TODO add more columns
attrs = [
"num_steps", "final_energy", "final_pressure",
#"final_min_force", "final_max_force",
#"ecut", "pawecutdg", "tsmear", "nkpt", "nsppol", "nspinor", "nspden",
"final_fmin", "final_fmax", "final_fmean", "final_fstd", "final_drift",
"initial_fmin", "initial_fmax", "initial_fmean", "initial_fstd", "initial_drift",
# TODO add more columns but must update HIST file
#"nsppol", "nspinor", "nspden",
#"ecut", "pawecutdg", "tsmear", "nkpt",
] + kwargs.pop("attrs", [])
rows, row_names = [], []
@ -439,15 +544,18 @@ class HistRobot(Robot):
row_names.append(label)
d = OrderedDict()
#fstas_dict = self.get_fstats_dict(step=-1)
initial_fstas_dict = hist.get_fstats_dict(step=0)
final_fstas_dict = hist.get_fstats_dict(step=-1)
# Add info on structure.
if with_geo:
d.update(hist.final_structure.get_dict4pandas(with_spglib=with_spglib))
for aname in attrs:
if aname in ("final_min_force", "final_max_force"):
value = fstas_dict[aname]
if aname in ("final_fmin", "final_fmax", "final_fmean", "final_fstd", "final_drift",):
value = final_fstas_dict[aname.replace("final_", "")]
elif aname in ("initial_fmin", "initial_fmax", "initial_fmean", "initial_fstd", "initial_drift"):
value = initial_fstas_dict[aname.replace("initial_", "")]
else:
value = getattr(hist, aname, None)
d[aname] = value
@ -461,6 +569,97 @@ class HistRobot(Robot):
index = row_names if index is None else index
return pd.DataFrame(rows, index=index, columns=list(rows[0].keys()))
@property
def what_list(self):
"""List with all quantities that can be plotted (what argument)."""
return ["energy", "abc", "angles", "volume", "pressure", "forces"]
@add_fig_kwargs
def gridplot(self, what="abc", sharex=False, sharey=False, fontsize=8, **kwargs):
"""
Plot the `what` value extracted from multiple HIST files on a grid.
Args:
what: Quantity to plot. Must be in ["energy", "abc", "angles", "volume", "pressure", "forces"]
sharex: True if xaxis should be shared.
sharey: True if yaxis should be shared.
fontsize: fontsize for legend.
Returns:
matplotlib figure.
"""
num_plots, ncols, nrows = len(self), 1, 1
if num_plots > 1:
ncols = 2
nrows = (num_plots // ncols) + (num_plots % ncols)
ax_list, fig, plt = get_axarray_fig_plt(None, nrows=nrows, ncols=ncols,
sharex=sharex, sharey=sharey, squeeze=False)
ax_list = ax_list.ravel()
for i, (ax, hist) in enumerate(zip(ax_list, self.abifiles)):
hist.plot_ax(ax, what, fontsize=fontsize, marker="o")
ax.set_title(hist.relpath, fontsize=fontsize)
ax.grid(True)
if i == len(ax_list) - 1:
ax.set_xlabel('Step')
else:
ax.set_xlabel('')
# Get around a bug in matplotlib.
if num_plots % ncols != 0:
ax_list[-1].plot([0, 1], [0, 1], lw=0)
ax_list[-1].axis('off')
return fig
@add_fig_kwargs
def combiplot(self, what_list=None, cmap="jet", fontsize=6, **kwargs):
"""
Plot multiple HIST files on a grid. One plot for each `what` value.
Args:
what_list: List of strings with the quantities to plot.
If None, all quanties are plotted.
cmap: matplotlib color map.
fontsize: fontisize for legend.
Returns:
matplotlib figure.
"""
what_list = (list_strings(what_list) if what_list is not None
else ["energy", "a", "b", "c", "alpha", "beta", "gamma", "volume", "pressure"])
num_plots, ncols, nrows = len(what_list), 1, 1
if num_plots > 1:
ncols = 2
nrows = (num_plots // ncols) + (num_plots % ncols)
ax_list, fig, plt = get_axarray_fig_plt(None, nrows=nrows, ncols=ncols,
sharex=True, sharey=False, squeeze=False)
ax_list = ax_list.ravel()
cmap = plt.get_cmap(cmap)
for i, (ax, what) in enumerate(zip(ax_list, what_list)):
for ih, hist in enumerate(self.abifiles):
label= None if i != 0 else hist.relpath
hist.plot_ax(ax, what, color=cmap(ih / len(self)), label=label, fontsize=fontsize)
if label is not None:
ax.legend(loc="best", fontsize=fontsize, shadow=True)
if i == len(ax_list) - 1:
ax.set_xlabel("Step")
else:
ax.set_xlabel("")
# Get around a bug in matplotlib.
if num_plots % ncols != 0:
ax_list[-1].plot([0, 1], [0, 1], lw=0)
ax_list[-1].axis('off')
return fig
def write_notebook(self, nbpath=None):
"""
Write a jupyter notebook to nbpath. If nbpath is None, a temporay file in the current
@ -472,7 +671,8 @@ class HistRobot(Robot):
nb.cells.extend([
#nbv.new_markdown_cell("# This is a markdown cell"),
nbv.new_code_cell("robot = abilab.HistRobot(*%s)\nrobot.trim_paths()\nrobot" % str(args)),
nbv.new_code_cell("df = robot.get_dataframe()\ndisplay(df)"),
nbv.new_code_cell("robot.get_dataframe()"),
nbv.new_code_cell("for what in robot.what_list: robot.gridplot(what=what, tight_layout=True);"),
])
# Mixins
@ -505,7 +705,7 @@ class HistReader(ETSF_Reader):
if num_pseudos != ntypat:
raise NotImplementedError("Alchemical mixing is not supported, num_pseudos != ntypat")
znucl, typat = self.read_value("znucl"), self.read_value("typat")
znucl, typat = self.read_value("znucl"), self.read_value("typat").astype(int)
#print(znucl.dtype, typat)
cart_forces_step = self.read_cart_forces(unit="eV ang^-1")

5
abipy/dynamics/tests/test_hist.py
View File

@ -93,5 +93,10 @@ class HistFileTest(AbipyTest):
df = robot.get_dataframe()
assert "angle1" in df
if self.has_matplotlib():
for what in robot.what_list:
assert robot.gridplot(what=what, show=False)
assert robot.combiplot(show=False)
if self.has_nbformat():
robot.write_notebook(nbpath=self.get_tmpname(text=True))

17
abipy/electrons/bse.py
View File

@ -571,7 +571,7 @@ class MdfPlotter(object):
self._mdfs[label] = mdf
@add_fig_kwargs
def plot(self, ax=None, cplx_mode="Im", qpoint=None, xlims=None, ylims=None, **kwargs):
def plot(self, ax=None, cplx_mode="Im", qpoint=None, xlims=None, ylims=None, fontsize=12, **kwargs):
"""
Get a matplotlib plot showing the MDFs.
@ -584,6 +584,7 @@ class MdfPlotter(object):
xlims: Set the data limits for the y-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
ylims: Same meaning as `ylims` but for the y-axis
fontsize: Legend and label fontsize.
"""
ax, fig, plt = get_ax_fig_plt(ax)
ax.grid(True)
@ -603,7 +604,7 @@ class MdfPlotter(object):
legends.append(r"%s: %s, %s $\varepsilon$" % (cmode, qtag, label))
# Set legends.
ax.legend(lines, legends, loc='best', shadow=False)
ax.legend(lines, legends, loc='best', fontsize=fontsize, shadow=True)
set_axlims(ax, xlims, "x")
set_axlims(ax, ylims, "y")
@ -802,8 +803,8 @@ class MultipleMdfPlotter(object):
# return fig
@add_fig_kwargs
def plot_mdftype_cplx(self, mdf_type, cplx_mode, qpoint=None, ax=None,
xlims=None, ylims=None, with_legend=True, with_xlabel=True, with_ylabel=True, **kwargs):
def plot_mdftype_cplx(self, mdf_type, cplx_mode, qpoint=None, ax=None, xlims=None, ylims=None,
with_legend=True, with_xlabel=True, with_ylabel=True, fontsize=12, **kwargs):
"""
Helper function to plot data corresponds to `mdf_type`, `cplx_mode`, `qpoint`.
@ -819,6 +820,7 @@ class MultipleMdfPlotter(object):
with_legend: True if legend should be added
with_xlabel:
with_ylabel:
fontsize: Legend and label fontsize.
Return: matplotlib figure
"""
@ -850,7 +852,7 @@ class MultipleMdfPlotter(object):
# Set legends.
if with_legend:
ax.legend(lines, legends, loc='best', shadow=False)
ax.legend(lines, legends, loc='best', fontsize=fontsize, shadow=True)
return fig
@ -910,11 +912,14 @@ class MdfRobot(Robot, RobotWithEbands):
"""
EXT = "MDF"
def plot(self, **kwargs):
"""Wraps plot method of `MultipleMdfPlotter`. kwargs passed to plot."""
return self.get_multimdf_plotter().plot(**kwargs)
def get_multimdf_plotter(self, cls=None):
"""
Return an instance of MultipleMdfPlotter to compare multiple dielectric functions.
"""
from abipy.electrons.bse import MultipleMdfPlotter
plotter = MultipleMdfPlotter() if cls is None else cls()
for label, mdf in self:

87
abipy/electrons/ebands.py
View File

@ -33,7 +33,6 @@ from abipy.iotools import ETSF_Reader
from abipy.tools import gaussian, duck
from abipy.tools.plotting import set_axlims, add_fig_kwargs, get_ax_fig_plt, get_ax3d_fig_plt
import logging
logger = logging.getLogger(__name__)
@ -624,7 +623,8 @@ class ElectronBands(Has_Structure):
with_spglib: If True, spglib is invoked to get the spacegroup symbol and number
"""
odict = OrderedDict([
("nsppol", self.nsppol), ("nkpt", self.nkpt), ("nband", self.nband_sk.min()),
("nsppol", self.nsppol), ("nspinor", self.nspinor), ("nspden", self.nspden),
("nkpt", self.nkpt), ("nband", self.nband_sk.min()),
("nelect", self.nelect), ("fermie", self.fermie),
])
@ -1455,7 +1455,7 @@ class ElectronBands(Has_Structure):
@add_fig_kwargs
def plot_ejdosvc(self, vrange, crange, method="gaussian", step=0.1, width=0.2, colormap="jet",
cumulative=True, ax=None, alpha=0.7, **kwargs):
cumulative=True, ax=None, alpha=0.7, fontsize=12, **kwargs):
"""
Plot the decomposition of the joint-density of States (JDOS).
@ -1473,6 +1473,7 @@ class ElectronBands(Has_Structure):
http://matplotlib.sourceforge.net/examples/pylab_examples/show_colormaps.html
cumulative: True for cumulative plots (default).
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: fontsize for legends and titles
Returns:
`matplotlib` figure
@ -1520,7 +1521,7 @@ class ElectronBands(Has_Structure):
tot_jdos.plot_ax(ax, color="k", lw=lw, label=r"Total JDOS, $\sigma=%s$" % s)
ax.legend(loc="best")
ax.legend(loc="best", shadow=True, fontsize=fontsize)
return fig
def apply_scissors(self, scissors):
@ -1820,8 +1821,16 @@ class ElectronBands(Has_Structure):
def to_xmgrace(self, filepath):
"""
Write xmgrace file with band structure energies and labels for high-symmetry k-points.
Args:
filepath: String with filename or stream.
"""
f = open(filepath, "wt")
is_stream = hasattr(filepath, "write")
if is_stream:
f = filepath
else:
f = open(filepath, "wt")
def w(s):
f.write(s)
f.write("\n")
@ -1885,7 +1894,8 @@ class ElectronBands(Has_Structure):
w('%d %.8E' % (ik, emef[spin, ik, band]))
w('&')
f.close()
if not is_stream:
f.close()
def to_bxsf(self, filepath):
"""
@ -2268,7 +2278,7 @@ class ElectronBandsPlotter(NotebookWriter):
return "\n\n".join(text)
@add_fig_kwargs
def combiplot(self, e0="fermie", ylims=None, width_ratios=(2, 1), **kwargs):
def combiplot(self, e0="fermie", ylims=None, width_ratios=(2, 1), fontsize=8, **kwargs):
"""
Plot the band structure and the DOS on the same figure.
Use `gridplot` to plot band structures on different figures.
@ -2289,6 +2299,7 @@ class ElectronBandsPlotter(NotebookWriter):
or scalar e.g. `left`. If left (right) is None, default values are used
width_ratios: Defines the ratio between the band structure plot and the dos plot.
Used when there are DOS stored in the plotter.
fontsize: fontsize for titles and legend.
Returns:
matplotlib figure.
@ -2341,7 +2352,7 @@ class ElectronBandsPlotter(NotebookWriter):
if i == 0:
ebands.decorate_ax(ax1)
ax1.legend(lines, legends, loc='upper right', shadow=True)
ax1.legend(lines, legends, loc='upper right', fontsize=fontsize, shadow=True)
# Add DOSes
if self.edoses_dict:
@ -2360,7 +2371,7 @@ class ElectronBandsPlotter(NotebookWriter):
return self.combiplot(*args, **kwargs)
@add_fig_kwargs
def gridplot(self, e0="fermie", with_dos=True, ylims=None, **kwargs):
def gridplot(self, e0="fermie", with_dos=True, ylims=None, fontsize=8, **kwargs):
"""
Plot multiple electron bandstructures and optionally DOSes on a grid.
@ -2386,6 +2397,7 @@ class ElectronBandsPlotter(NotebookWriter):
with_dos: True if DOS should be printed.
ylims: Set the data limits for the y-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
fontsize: fontsize for titles and legend
Returns:
matplotlib figure.
@ -2410,7 +2422,7 @@ class ElectronBandsPlotter(NotebookWriter):
for i, (ebands, ax) in enumerate(zip(ebands_list, axes)):
ebands.plot(ax=ax, e0=e0, show=False)
set_axlims(ax, ylims, "y")
if titles is not None: ax.set_title(titles[i])
if titles is not None: ax.set_title(titles[i], fontsize=fontsize)
if i % ncols != 0:
ax.set_ylabel("")
@ -2432,7 +2444,7 @@ class ElectronBandsPlotter(NotebookWriter):
mye0 = ebands.get_e0(e0) if e0 != "edos_fermie" else edos.fermie
ebands.plot_with_edos(edos, e0=mye0, axlist=(ax1, ax2), show=False)
if titles is not None: ax1.set_title(titles[i])
if titles is not None: ax1.set_title(titles[i], fontsize=fontsize)
if i % ncols != 0:
for ax in (ax1, ax2):
ax.set_ylabel("")
@ -2440,7 +2452,7 @@ class ElectronBandsPlotter(NotebookWriter):
return fig
@add_fig_kwargs
def boxplot(self, e0="fermie", brange=None, swarm=False, **kwargs):
def boxplot(self, e0="fermie", brange=None, swarm=False, fontsize=8, **kwargs):
"""
Use seaborn to draw a box plot to show distributions of eigenvalues with respect to the band index.
Band structures are drawn on different subplots.
@ -2452,6 +2464,7 @@ class ElectronBandsPlotter(NotebookWriter):
- None: Don't shift energies, equivalent to e0=0
brange: Only bands such as `brange[0] <= band_index < brange[1]` are included in the plot.
swarm: True to show the datapoints on top of the boxes
fontsize: Fontsize for title.
kwargs: Keyword arguments passed to seaborn boxplot.
"""
# Build grid of plots.
@ -2468,7 +2481,7 @@ class ElectronBandsPlotter(NotebookWriter):
for (label, ebands), ax in zip(self.ebands_dict.items(), ax_list):
ebands.boxplot(ax=ax, brange=brange, show=False)
ax.set_title(label)
ax.set_title(label, fontsize=fontsize)
return fig
@ -3102,7 +3115,7 @@ class ElectronDosPlotter(NotebookWriter):
self.edoses_dict[label] = ElectronDos.as_edos(edos, edos_kwargs)
@add_fig_kwargs
def combiplot(self, ax=None, e0="fermie", xlims=None, **kwargs):
def combiplot(self, ax=None, e0="fermie", xlims=None, fontsize=8, **kwargs):
"""
Plot the the DOSes on the same figure.
Use `gridplot` to plot DOSes on different figures.
@ -3115,6 +3128,7 @@ class ElectronDosPlotter(NotebookWriter):
- None: Don't shift energies, equivalent to e0=0
xlims: Set the data limits for the x-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
fontsize: fontsize for titles and legend
Returns:
`matplotlib` figure.
@ -3134,7 +3148,7 @@ class ElectronDosPlotter(NotebookWriter):
ax.set_xlabel("Energy [eV]")
ax.set_ylabel('DOS [states/eV]')
set_axlims(ax, xlims, "x")
ax.legend(loc="best")
ax.legend(loc="best", shadow=True, fontsize=fontsize)
return fig
@ -3435,13 +3449,18 @@ class Bands3D(object):
Require k-points in IBZ and gamma-centered k-mesh.
Args:
filepath: BXSF filename.
filepath: BXSF filename or stream.
unit: Input energies are in unit `unit`.
"""
from abipy.iotools import bxsf_write
with open(filepath, "wt") as fh:
bxsf_write(fh, self.structure, self.nsppol, self.nband, self.kdivs, self.ucdata_sbk, self.fermie, unit=unit)
return filepath
if hasattr(filepath, "write"):
return bxsf_write(filepath, self.structure, self.nsppol, self.nband, self.kdivs,
self.ucdata_sbk, self.fermie, unit=unit)
else:
with open(filepath, "wt") as fh:
bxsf_write(fh, self.structure, self.nsppol, self.nband, self.kdivs,
self.ucdata_sbk, self.fermie, unit=unit)
return filepath
def get_e0(self, e0):
"""
@ -3567,7 +3586,7 @@ class Bands3D(object):
return figure
@add_fig_kwargs
def plot_contour(self, band, spin=0, plane="xy", elevation=0, ax=None, **kwargs):
def plot_contour(self, band, spin=0, plane="xy", elevation=0, ax=None, fontsize=10, **kwargs):
"""
Contour plot with matplotlib.
@ -3577,6 +3596,7 @@ class Bands3D(object):
plane:
elevation:
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: Label and title fontsize.
Return:
`matplotlib` figure
@ -3590,9 +3610,9 @@ class Bands3D(object):
ax, fig, plt = get_ax_fig_plt(ax=ax)
x, y = np.meshgrid(x, y)
c = ax.contour(x, y, fxy, **kwargs)
ax.clabel(c, inline=1, fontsize=10)
ax.clabel(c, inline=1, fontsize=fontsize)
kvert = dict(xy="z", xz="y", yz="x")[plane]
ax.set_title("Band %s in %s plane at $K_%s=%d$" % (band, plane, kvert, elevation))
ax.set_title("Band %s in %s plane at $K_%s=%d$" % (band, plane, kvert, elevation), fontsize=fontsize)
ax.grid(True)
ax.set_xlabel("$K_%s$" % plane[0])
ax.set_ylabel("$K_%s$" % plane[1])
@ -3668,6 +3688,29 @@ class RobotWithEbands(object):
"""
Mixin class for robots associated to files with `ElectronBands`.
"""
def combiplot_ebands(self, **kwargs):
"""Wraps combiplot method of `ElectronBandsPlotter`. kwargs passed to combiplot."""
return self.get_ebands_plotter().combiplot(**kwargs)
def gridplot_ebands(self, **kwargs):
"""Wraps gridplot method of `ElectronBandsPlotter`. kwargs passed to gridplot."""
return self.get_ebands_plotter().gridplot(**kwargs)
def boxplot_ebands(self, **kwargs):
"""Wraps boxplot method of `ElectronBandsPlotter`. kwargs passed to boxplot."""
return self.get_ebands_plotter().boxplot(**kwargs)
def combiboxplot_ebands(self, **kwargs):
"""Wraps combiboxplot method of `ElectronDosPlotter`. kwargs passed to combiboxplot."""
return self.get_ebands_plotter().combiboxplot(**kwargs)
def combiplot_edos(self, **kwargs):
"""Wraps combiplot method of `ElectronDosPlotter`. kwargs passed to combiplot."""
return self.get_edos_plotter().combiplot(**kwargs)
def gridplot_edos(self, **kwargs):
"""Wraps gridplot method of `ElectronDosPlotter`. kwargs passed to gridplot."""
return self.get_edos_plotter().gridplot(**kwargs)
def get_ebands_plotter(self, filter_abifile=None, cls=None):
"""

16
abipy/electrons/fatbands.py
View File

@ -938,7 +938,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
@add_fig_kwargs
def plot_pjdos_lview(self, e0="fermie", lmax=None, method="gaussian", step=0.1, width=0.2,
stacked=True, combined_spins=True, axmat=None, exchange_xy=False,
with_info=True, with_spin_sign=True, xlims=None, ylims=None, **kwargs):
with_info=True, with_spin_sign=True, xlims=None, ylims=None, fontsize=12, **kwargs):
"""
Plot the PJ-DOS on a linear mesh.
@ -961,6 +961,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
xlims: Set the data limits for the x-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
ylims: Same meaning as `xlims` but for the y-axis
fontsize: Legend and label fontsize
Returns:
`matplotlib` figure
@ -1067,7 +1068,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
else:
title = "%s, %s" % (self.l2tex[l], self.spin2tex[spin]) if self.nsppol == 2 else \
self.l2tex[l]
ax.set_title(title)
ax.set_title(title, fontsize=fontsize)
ax.grid(True)
# Display yticklabels on the first plot and last plot only.
@ -1075,7 +1076,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
ax.set_xlabel("Energy [eV]")
if l == 0:
if with_info:
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
if exchange_xy:
ax.set_xlabel('DOS [states/eV]')
else:
@ -1096,7 +1097,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
@add_fig_kwargs
def plot_pjdos_typeview(self, e0="fermie", lmax=None, method="gaussian", step=0.1, width=0.2,
stacked=True, combined_spins=True, axmat=None, exchange_xy=False,
with_info=True, with_spin_sign=True, xlims=None, ylims=None, **kwargs):
with_info=True, with_spin_sign=True, xlims=None, ylims=None, fontsize=12, **kwargs):
"""
Plot the PJ-DOS on a linear mesh.
@ -1119,6 +1120,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
xlims: Set the data limits for the x-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
ylims: Same meaning as `xlims` but for the y-axis
fontsize: Legend and label fontsize.
Returns:
`matplotlib` figure
@ -1224,7 +1226,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
title = "Type: %s" % symbol
else:
title = "%s, %s" % (symbol, self.spin2tex[spin]) if self.nsppol == 2 else symbol
ax.set_title(title)
ax.set_title(title, fontsize=fontsize)
ax.grid(True)
# Display yticklabels on the first plot and last plot only.
@ -1232,7 +1234,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
ax.set_xlabel("Energy [eV]")
if itype == 0:
if with_info:
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
if exchange_xy:
ax.set_xlabel('DOS [states/eV]')
else:
@ -1529,7 +1531,7 @@ class FatBandsFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, N
# ax.grid(True)
# ax.set_xlabel("Energy [eV]")
# set_axlims(ax, xlims, "x")
# ax.legend(loc="best")
# ax.legend(loc="best", fontsize=fontsize, shadow=True)
# return fig

5
abipy/electrons/fold2bloch.py
View File

@ -171,7 +171,7 @@ class Fold2BlochNcfile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBand
@add_fig_kwargs
def plot_unfolded(self, kbounds, klabels, ylims=None, dist_tol=1e-12, verbose=0,
colormap="afmhot", facecolor="black", ax=None, **kwargs):
colormap="afmhot", facecolor="black", ax=None, fontsize=12, **kwargs):
r"""
Plot unfolded band structure with spectral weights.
@ -187,6 +187,7 @@ class Fold2BlochNcfile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBand
http://matplotlib.sourceforge.net/examples/pylab_examples/show_colormaps.html
facecolor:
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: Legend and title fontsize.
Returns:
`matplotlib` figure
@ -224,7 +225,7 @@ class Fold2BlochNcfile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBand
ax.grid(True)
ax.set_ylabel('Energy [eV]')
set_axlims(ax, ylims, "y")
if self.nss == 2: ax.legend(loc="best")
if self.nss == 2: ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig

158
abipy/electrons/gsr.py
View File

@ -9,13 +9,14 @@ import pymatgen.core.units as units
from collections import OrderedDict, Iterable, defaultdict
from tabulate import tabulate
from monty.string import is_string, list_strings, marquee
from monty.collections import AttrDict
from monty.termcolor import cprint
from monty.collections import AttrDict, dict2namedtuple
from monty.functools import lazy_property
from pymatgen.core.units import EnergyArray, ArrayWithUnit
from pymatgen.entries.computed_entries import ComputedEntry, ComputedStructureEntry
from abipy.core.mixins import AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, NotebookWriter
from prettytable import PrettyTable
from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt
from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, get_axarray_fig_plt
from abipy.abio.robots import Robot
from abipy.electrons.ebands import ElectronsReader, RobotWithEbands
@ -118,12 +119,12 @@ class GsrFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, Notebo
@lazy_property
def energy(self):
"""Total energy"""
"""Total energy in eV"""
return units.Energy(self.reader.read_value("etotal"), "Ha").to("eV")
@lazy_property
def energy_per_atom(self):
"""Total energy / number_of_atoms"""
"""Total energy / number_of_atoms (eV units)"""
return self.energy / len(self.structure)
@lazy_property
@ -397,6 +398,11 @@ class GsrRobot(Robot, RobotWithEbands):
for label, gsr in self:
row_names.append(label)
d = OrderedDict()
# Add info on structure.
if with_geo:
d.update(gsr.structure.get_dict4pandas(with_spglib=True))
for aname in attrs:
if aname == "nkpt":
value = len(gsr.ebands.kpoints)
@ -405,10 +411,6 @@ class GsrRobot(Robot, RobotWithEbands):
if value is None: value = getattr(gsr.ebands, aname, None)
d[aname] = value
# Add info on structure.
if with_geo:
d.update(gsr.structure.get_dict4pandas(with_spglib=True))
# Execute functions
if funcs is not None: d.update(self._exec_funcs(funcs, gsr))
rows.append(d)
@ -416,42 +418,92 @@ class GsrRobot(Robot, RobotWithEbands):
row_names = row_names if not abspath else self._to_relpaths(row_names)
return pd.DataFrame(rows, index=row_names, columns=list(rows[0].keys()))
# FIXME: EOS has been changed in pymatgen.
def eos_fit(self, eos_name="murnaghan"):
def get_eos_fits_dataframe(self, eos_names="murnaghan"):
"""
Fit energy as function of volume to get the equation of state, equilibrium volume,
bulk modulus and its derivative wrt to pressure.
Fit energy as function of volume to get the equation of state,
equilibrium volume, bulk modulus and its derivative wrt to pressure.
Args:
eos_name:
eos_names: String or list of strings with EOS names.
For the list of available models, see pymatgen.analysis.eos.
Return:
(fits, dataframe) namedtuple.
fits is a list of `EOSFit object`
dataframe is a pandas Dataframe with the final results.
"""
# Read volumes and energies from the GSR files.
energies, volumes = [], []
for label, gsr in self:
energies.append(gsr.energy)
volumes.append(gsr.structure.volume)
energies.append(float(gsr.energy))
volumes.append(float(gsr.structure.volume))
# Order data by volumes if needed.
if np.any(np.diff(volumes) < 0):
ves = sorted(zip(volumes, energies), key=lambda t: t[0])
volumes = [t[0] for t in ves]
energies = [t[1] for t in ves]
# Note that eos.fit expects lengths in Angstrom, and energies in eV.
# I'm also monkey-patching the plot method.
if eos_name != "all":
fit = EOS(eos_name=eos_name).fit(volumes, energies)
fit.plot = _my_fit_plot
return fit
else:
from pymatgen.analysis.eos import EOS
if eos_names == "all":
# Use all the available models.
fits, rows = [], []
models = list(EOS.MODELS.keys())
for eos_name in models:
fit = EOS(eos_name=eos_name).fit(volumes, energies)
fit.plot = _my_fit_plot
fits.append(fit)
rows.append(OrderedDict([(aname, getattr(fit, aname)) for aname in
("v0", "e0", "b0_GPa", "b1")]))
eos_names = [n for n in EOS.MODELS if n not in ("deltafactor", "numerical_eos")]
else:
eos_names = list_strings(eos_names)
frame = pd.DataFrame(rows, index=models, columns=list(rows[0].keys()))
return fits, frame
#return dict2namedtuple(fits=fits, frame=frame)
fits, index, rows = [], [], []
for eos_name in eos_names:
try:
fit = EOS(eos_name=eos_name).fit(volumes, energies)
except Exception as exc:
cprint("EOS %s raised exception:\n%s" % (eos_name, str(exc)))
continue
# Replace plot with plot_ax method
fit.plot = fit.plot_ax
fits.append(fit)
index.append(eos_name)
rows.append(OrderedDict([(aname, getattr(fit, aname)) for aname in
("v0", "e0", "b0_GPa", "b1")]))
dataframe = pd.DataFrame(rows, index=index, columns=list(rows[0].keys()) if rows else None)
return dict2namedtuple(fits=fits, dataframe=dataframe)
@add_fig_kwargs
def gridplot_eos(self, eos_names="all", fontsize=6, **kwargs):
"""
Plot multiple EOS on a grid with captions showing the final results.
Args:
eos_names: String or list of strings with EOS names. See pymatgen.analysis.EOS
fontsize: Fontsize used for caption text.
Returns:
matplotlib figure.
"""
r = self.get_eos_fits_dataframe(eos_names=eos_names)
num_plots, ncols, nrows = len(r.fits), 1, 1
if num_plots > 1:
ncols = 2
nrows = (num_plots // ncols) + (num_plots % ncols)
# Build grid of plots.
ax_list, fig, plt = get_axarray_fig_plt(None, nrows=nrows, ncols=ncols,
sharex=False, sharey=False, squeeze=False)
ax_list = ax_list.ravel()
for i, (fit, ax) in enumerate(zip(r.fits, ax_list)):
fit.plot_ax(ax=ax, fontsize=fontsize, label="", show=False)
# Get around a bug in matplotlib
if num_plots % ncols != 0:
ax_list[-1].plot([0, 1], [0, 1], lw=0)
ax_list[-1].axis('off')
return fig
#def get_phasediagram_results(self):
# from abipy.core.restapi import PhaseDiagramResults
@ -477,6 +529,7 @@ class GsrRobot(Robot, RobotWithEbands):
nbv.new_code_cell("#anim = ebands_plotter.animate();"),
nbv.new_code_cell("edos_plotter = robot.get_edos_plotter()"),
nbv.new_code_cell("edos_plotter.ipw_select_plot()"),
nbv.new_code_cell("#robot.gridplot_eos();"),
])
# Mixins
@ -484,46 +537,3 @@ class GsrRobot(Robot, RobotWithEbands):
nb.cells.extend(self.get_ebands_code_cells())
return self._write_nb_nbpath(nb, nbpath)
@add_fig_kwargs
def _my_fit_plot(self, ax=None, **kwargs):
"""
Plot the equation of state.
Args:
Returns:
Matplotlib figure object.
"""
ax, fig, plt = get_ax_fig_plt(ax=ax)
color = kwargs.get("color", "r")
label = kwargs.get("label", "{} fit".format(self.__class__.__name__))
lines = ["Equation of State: %s" % self.__class__.__name__,
"Minimum energy = %1.2f eV" % self.e0,
"Minimum or reference volume = %1.2f Ang^3" % self.v0,
"Bulk modulus = %1.2f eV/Ang^3 = %1.2f GPa" %
(self.b0, self.b0_GPa),
"Derivative of bulk modulus wrt pressure = %1.2f" % self.b1]
text = "\n".join(lines)
text = kwargs.get("text", text)
# Plot input data.
ax.plot(self.volumes, self.energies, linestyle="None", marker="o", color=color)
# Plot eos fit.
vmin, vmax = min(self.volumes), max(self.volumes)
vmin, vmax = (vmin - 0.01 * abs(vmin), vmax + 0.01 * abs(vmax))
vfit = np.linspace(vmin, vmax, 100)
ax.plot(vfit, self.func(vfit), linestyle="dashed", color=color, label=label)
ax.grid(True)
ax.xlabel("Volume $\\AA^3$")
ax.ylabel("Energy (eV)")
ax.legend(loc="best", shadow=True)
# Add text with fit parameters.
ax.text(0.4, 0.5, text, transform=ax.transAxes)
return fig

173
abipy/electrons/gw.py
View File

@ -8,14 +8,13 @@ import numpy as np
import pandas as pd
from collections import namedtuple, OrderedDict, Iterable, defaultdict
from six.moves import cStringIO
from monty.string import list_strings, is_string, marquee
from monty.collections import AttrDict, dict2namedtuple
from monty.functools import lazy_property
from monty.termcolor import cprint
from monty.dev import deprecated
from monty.bisect import find_le, find_ge
from prettytable import PrettyTable
from six.moves import cStringIO
from abipy.core.func1d import Function1D
from abipy.core.kpoints import Kpoint, KpointList, Kpath, IrredZone, has_timrev_from_kptopt
from abipy.core.mixins import AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter
@ -33,6 +32,7 @@ __all__ = [
"QPState",
"SigresFile",
"SigresPlotter",
"SigresRobot",
]
@ -239,10 +239,12 @@ class QPList(list):
"""Return an arrays with the :class:`QPState` corrections."""
return self.get_field("qpeme0")
@deprecated(message="to_table is deprecated and will be removed in v0.4")
def to_table(self):
"""Return a table (list of list of strings)."""
header = QPState.get_fields(exclude=["spin", "kpoint"])
# TODO: Use tabulate or pd
from prettytable import PrettyTable
table = PrettyTable(header)
for qp in self:
@ -438,8 +440,8 @@ class Sigmaw(object):
self.xc = Function1D(self.wmesh, sigmaxc_values)
self.spfunc = Function1D(self.wmesh, spfunc_values)
def plot_ax(self, ax, w="a", **kwargs):
"""Helper function to plot data on the axis ax."""
def plot_ax(self, ax, w="a", fontsize=12, **kwargs):
"""Helper function to plot data on the axis ax with fontsize"""
#if not kwargs: kwargs = {"color": "black", "linewidth": 2.0}
lines = []
@ -450,7 +452,7 @@ class Sigmaw(object):
label = kwargs.get("label", r"$\Sigma(\omega)$")
extend(f.plot_ax(ax, cplx_mode="re", label="Re " + label))
extend(f.plot_ax(ax, cplx_mode="im", label="Im " + label))
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
#ax.set_ylabel('Energy [eV]')
elif w == "a":
@ -461,7 +463,7 @@ class Sigmaw(object):
#ax2 = ax.twinx()
#extend(f.cumintegral().plot_ax(ax2, label="$I(\omega) = \int_{-\infty}^{\omega} A(\omega')d\omega'$"))
#ax.set_ylabel('Energy [eV]')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
else:
raise ValueError("Don't know how to handle what option %s" % w)
@ -685,12 +687,12 @@ class SigresPlotter(Iterable):
assert len(xvalues) == len(self)
self._xvalues = xvalues
def decorate_ax(self, ax, **kwargs):
def decorate_ax(self, ax, fontsize=12, **kwargs):
ax.grid(True)
if self.param_name is not None:
ax.set_xlabel(self.param_name)
ax.set_ylabel('Energy [eV]')
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
title = kwargs.pop("title", None)
if title is not None: ax.set_title(title)
@ -935,18 +937,33 @@ class SigresFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter)
return self.to_string()
def to_string(self, verbose=0):
"""String representation."""
"""String representation with verbosity level `verbose`."""
lines = []; app = lines.append
app(marquee("File Info", mark="="))
app(self.filestat(as_string=True))
app("")
app(self.ebands.to_string(title="Kohn-Sham bands"))
app(self.structure.to_string(verbose=verbose, title="Structure"))
app("")
app(self.ebands.to_string(title="Kohn-Sham bands", with_structure=False))
# TODO:
# Finalize the implementation: add GW quantities.
# Fix header.
# TODO: Finalize the implementation: add GW metadata.
app(marquee("QP direct gaps in eV", mark="="))
for kgw in self.gwkpoints:
for spin in range(self.nsppol):
qp_dirgap = self.get_qpgap(spin, kgw)
#ks_dirgap =
app("QP_dirgap: %.3f for K-point: %s, spin: %s" % (qp_dirgap, repr(kgw), spin))
strio = cStringIO()
self.print_qps(file=strio)
strio.seek(0)
app("")
app(marquee("QP results for each k-point and spin (All in eV)", mark="="))
app("".join(strio))
app("")
# TODO: Fix header.
#if verbose > 1:
# app("")
# app(marquee("Abinit Header", mark="="))
@ -1005,9 +1022,22 @@ class SigresFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter)
wmesh, spf_values = self.reader.read_spfunc(spin, kpoint, band)
return Function1D(wmesh, spf_values)
@deprecated(message="print_qps is deprecated and will be removed in version 0.4")
def print_qps(self, **kwargs):
self.reader.print_qps(**kwargs)
def print_qps(self, precision=3, ignore_imag=True, file=sys.stdout):
"""
Print QP results to stream `file`.
Args:
precision: Number of significant digits.
ignore_imag: True if imaginary part should be ignored.
file: Output stream.
"""
from abipy.tools.printing import print_dataframe
keys = "band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0".split()
for gwkpoint in self.gwkpoints:
for spin in range(self.nsppol):
df_sk = self.get_dataframe_sk(spin, gwkpoint, ignore_imag=ignore_imag)[keys]
print_dataframe(df_sk, title="K-point: %s, spin: %s" % (repr(gwkpoint), spin),
precision=precision, file=file)
@add_fig_kwargs
def plot_qps_vs_e0(self, with_fields="all", exclude_fields=None, axlist=None, label=None, **kwargs):
@ -1135,22 +1165,25 @@ class SigresFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter)
# FIXME: To maintain previous interface.
to_dataframe = get_dataframe
def get_dataframe_sk(self, spin, kpoint, index=None, ignore_imag=False):
def get_dataframe_sk(self, spin, kpoint, index=None, ignore_imag=False, with_params=True):
"""
Returns pandas DataFrame with QP results for the given (spin, k-point).
Args:
ignore_imag: Only real part is returned if `ignore_imag`.
with_params: True to include convergence paramenters.
"""
rows, bands = [], []
# FIXME start and stop should depend on k
for band in range(self.min_gwbstart, self.max_gwbstop):
# bstart and bstop depends on kpoint.
ik_gw = self.reader.gwkpt2seqindex(kpoint)
for band in range(self.gwbstart_sk[spin, ik_gw], self.gwbstop_sk[spin, ik_gw]):
bands.append(band)
# Build dictionary with the QP results.
qpstate = self.reader.read_qp(spin, kpoint, band, ignore_imag=ignore_imag)
d = qpstate.as_dict()
# Add other entries that may be useful when comparing different calculations.
d.update(self.params)
if with_params:
d.update(self.params)
rows.append(d)
import pandas as pd
@ -1326,7 +1359,7 @@ class SigresFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter)
eigens_kmesh = qp_corrs if only_corrections else ref_eigens + qp_corrs
# Build new ebands object with k-mesh
#ksampling = KSamplingInfo.from_mpdivs(mpdivs=kmesh, shifts=[0,0,0], kptopt=1)
#ksampling = KSamplingInfo.from_mpdivs(mpdivs=kmesh, shifts=[0, 0, 0], kptopt=1)
kpts_kmesh = IrredZone(self.structure.reciprocal_lattice, dos_kcoords, weights=dos_weights,
names=None, ksampling=ks_ebands_kmesh.kpoints.ksampling)
occfacts_kmesh = np.zeros(eigens_kmesh.shape)
@ -1633,7 +1666,7 @@ class SigresReader(ETSF_Reader):
qps = []
for gwkpoint in self.gwkpoints:
ik = self.gwkpt2seqindex(gwkpoint)
for band in range(self.gwbstart_sk[spin,ik], self.gwbstop_sk[spin,ik]):
for band in range(self.gwbstart_sk[spin, ik], self.gwbstop_sk[spin, ik]):
qps.append(self.read_qp(spin, gwkpoint, band, ignore_imag=ignore_imag))
qps_spin[spin] = QPList(qps)
@ -1664,7 +1697,8 @@ class SigresReader(ETSF_Reader):
Only real part is returned if `ignore_imag`.
"""
ik_file = self.kpt2fileindex(kpoint)
ib_file = band - self.gwbstart_sk[spin, self.gwkpt2seqindex(kpoint)]
# Must shift band index (see fortran code that allocates with mdbgw)
ib_gw = band - self.min_gwbstart
def ri(a):
return np.real(a) if ignore_imag else a
@ -1676,11 +1710,12 @@ class SigresReader(ETSF_Reader):
e0=self.read_e0(spin, ik_file, band),
qpe=ri(self._egw[spin, ik_file, band]),
qpe_diago=ri(self._en_qp_diago[spin, ik_file, band]),
vxcme=self._vxcme[spin, ik_file, ib_file],
sigxme=self._sigxme[spin, ik_file, ib_file],
sigcmee0=ri(self._sigcmee0[spin, ik_file, ib_file]),
vUme=self._vUme[spin, ik_file, ib_file],
ze0=ri(self._ze0[spin, ik_file, ib_file]),
# Note ib_gw index.
vxcme=self._vxcme[spin, ik_file, ib_gw],
sigxme=self._sigxme[spin, ik_file, ib_gw],
sigcmee0=ri(self._sigcmee0[spin, ik_file, ib_gw]),
vUme=self._vUme[spin, ik_file, ib_gw],
ze0=ri(self._ze0[spin, ik_file, ib_gw]),
)
def read_qpgaps(self):
@ -1710,14 +1745,16 @@ class SigresReader(ETSF_Reader):
raise ValueError("%s does not contain spectral function data" % self.path)
ik = self.kpt2fileindex(kpoint)
ib = band - self.gwbstart_sk[spin, self.gwkpt2seqindex(kpoint)]
# Must shift band index (see fortran code that allocates with mdbgw)
ib_gw = band - self.min_gwbstart
#ib_gw = band - self.gwbstart_sk[spin, self.gwkpt2seqindex(kpoint)]
aim_sigc = np.abs(self._sigcme[spin,:,ik,ib].imag)
aim_sigc = np.abs(self._sigcme[spin,:,ik,ib_gw].imag)
den = np.zeros(self.nomega_r)
for io, omega in enumerate(self._omega_r):
den[io] = (omega - self._hhartree[spin,ik,ib,ib].real - self._sigxcme[spin,io,ik,ib].real) ** 2 + \
self._sigcme[spin,io,ik,ib].imag ** 2
den[io] = (omega - self._hhartree[spin,ik,ib_gw,ib_gw].real - self._sigxcme[spin,io,ik,ib_gw].real) ** 2 + \
self._sigcme[spin,io,ik,ib_gw].imag ** 2
return self._omega_r, 1./np.pi * (aim_sigc/den)
@ -1741,57 +1778,18 @@ class SigresReader(ETSF_Reader):
"gwcalctyp", "scissor_ene",
]
# Read data and convert to scalar to avoid problems with pandas dataframes.
# Old sigres files may not have all the metadata.
params = AttrDict()
for pname in param_names:
params[pname] = self.read_value(pname, default=None)
v = self.read_value(pname, default=None)
params[pname] = v if v is None else np.asscalar(v)
# Other quantities that might be subject to convergence studies.
params["nkibz"] = len(self.ibz)
return params
@deprecated(message="print_qps is deprecated and will be removed in version 0.4")
def print_qps(self, spin=None, kpoints=None, bands=None, fmt=None, stream=sys.stdout):
"""
Args:
spin: Spin index, if None all spins are considered
kpoints: List of k-points to select. Default: all kpoints
bands: List of bands to select. Default is all bands
fmt: Format string passe to `to_strdict`
stream: file-like object.
Returns
List of tables.
"""
spins = range(self.nsppol) if spin is None else [spin]
kpoints = self.gwkpoints if kpoints is None else [kpoints]
if bands is not None: bands = [bands]
header = QPState.get_fields(exclude=["spin", "kpoint"])
tables = []
for spin in spins:
for kpoint in kpoints:
# TODO: Use tabulate or pd
table_sk = PrettyTable(header)
if bands is None:
ik = self.gwkpt2seqindex(kpoint)
bands = range(self.gwbstart_sk[spin,ik], self.gwbstop_sk[spin,ik])
for band in bands:
qp = self.read_qp(spin, kpoint, band)
d = qp.to_strdict(fmt=fmt)
table_sk.add_row([d[k] for k in header])
stream.write("\nkpoint: %s, spin: %s, energy units: eV (NB: bands start from zero)\n" % (kpoint, spin))
print(table_sk, file=stream)
stream.write("\n")
# Add it to tables.
tables.append(table_sk)
return tables
#def read_mel(self, mel_name, spin, kpoint, band, band2=None):
# array = self.read_value(mel_name)
@ -1825,7 +1823,7 @@ class SigresRobot(Robot, RobotWithEbands):
def get_qpgaps_dataframe(self, spin=None, kpoint=None, with_geo=False, abspath=False, funcs=None, **kwargs):
"""
Return a pandas DataFrame with the most important results for the given (spin, kpoint).
Return a pandas DataFrame with the QP gaps for all files in the robot.
Args:
spin: Spin index.
@ -1871,21 +1869,14 @@ class SigresRobot(Robot, RobotWithEbands):
row_names = row_names if not abspath else self._to_relpaths(row_names)
return pd.DataFrame(rows, index=row_names, columns=list(rows[0].keys()))
def plot_conv_qpgap(self, x_vars, show=True, **kwargs):
"""
Plot the convergence of the Quasi-particle gap.
kwargs are passed to :class:`seaborn.PairGrid`.
"""
import matplotlib.pyplot as plt
import seaborn.apionly as sns
# An alias to have a common API for robots.
get_dataframe = get_qpgaps_dataframe
data = self.get_qpgaps_dataframe()
#print(list(data.keys()))
grid = sns.PairGrid(data, x_vars=x_vars, y_vars="qpgap", **kwargs)
grid.map(plt.plot, marker="o")
grid.add_legend()
if show:
plt.show()
#@add_fig_kwargs
#def plot_qpgaps(self, ax=None, spin=None, kpoint=None, **kwargs):
#@add_fig_kwargs
#def plot_qpenes(self, spin=None, kpoint=None, band=None, hspan=0.01, **kwargs):
def write_notebook(self, nbpath=None):
"""
@ -1898,7 +1889,7 @@ class SigresRobot(Robot, RobotWithEbands):
nb.cells.extend([
#nbv.new_markdown_cell("# This is a markdown cell"),
nbv.new_code_cell("robot = abilab.SigresRobot(*%s)\nrobot.trim_paths()\nrobot" % str(args)),
#nbv.new_code_cell("df = robot.get_qpgaps_dataframe(spin=None, kpoint=None, with_geo=False, **kwargs)"),
#nbv.new_code_cell("robot.get_qpgaps_dataframe(spin=None, kpoint=None, with_geo=False, **kwargs)"),
#nbv.new_code_cell("plotter = robot.get_ebands_plotter()"),
])

12
abipy/electrons/optic.py
View File

@ -252,7 +252,7 @@ class OpticNcFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, No
@add_fig_kwargs
def plot_linear_epsilon(self, components="all", what="im", itemp=0,
ax=None, xlims=None, with_xlabel=True, label=None, **kwargs):
ax=None, xlims=None, with_xlabel=True, label=None, fontsize=12, **kwargs):
"""
Plot components of the linear dielectric function.
@ -266,6 +266,7 @@ class OpticNcFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, No
or scalar e.g. `left`. If left (right) is None, default values are used.
with_xlabel: True if x-label should be added.
label: True to add legend label to each curve.
fontsize: Legend and label fontsize.
Returns:
`matplotlib` figure
@ -283,7 +284,7 @@ class OpticNcFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, No
ax.grid(True)
if with_xlabel: ax.set_xlabel('Photon Energy [eV]')
set_axlims(ax, xlims, "x")
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@ -318,7 +319,7 @@ class OpticNcFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, No
@add_fig_kwargs
def plot_chi2(self, key, components="all", what="abs", itemp=0, decompose=False,
ax=None, xlims=None, with_xlabel=True, label=None, **kwargs):
ax=None, xlims=None, with_xlabel=True, label=None, fontsize=12, **kwargs):
"""
Low-level function to plot chi2 tensor.
@ -334,6 +335,7 @@ class OpticNcFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, No
or scalar e.g. `left`. If left (right) is None, default values are used.
with_xlabel: True to add x-label.
label: True to add legend label to each curve.
fontsize: Legend and label fontsize.
Returns:
`matplotlib` figure
@ -353,7 +355,7 @@ class OpticNcFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, No
ax.grid(True)
if with_xlabel: ax.set_xlabel('Photon Energy [eV]')
set_axlims(ax, xlims, "x")
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@ -488,7 +490,7 @@ class OpticRobot(Robot, RobotWithEbands):
available in each file. Use keys from ALL_CHIS.
"""
od = OrderedDict()
for ncfile in self.ncfiles:
for ncfile in self.abifiles:
for chiname in ALL_CHIS:
comps = ncfile.reader.computed_components[chiname]
if chiname not in od:

17
abipy/electrons/scr.py
View File

@ -152,7 +152,7 @@ class ScrFile(AbinitNcFile, Has_Header, Has_Structure, NotebookWriter):
return self.reader.read_params()
@add_fig_kwargs
def plot_emacro(self, cplx_mode="re-im", ax=None, xlims=None, **kwargs):
def plot_emacro(self, cplx_mode="re-im", ax=None, xlims=None, fontsize=12, **kwargs):
r"""
Plot the macroscopic dielectric function with local-field effects.
@ -165,6 +165,7 @@ class ScrFile(AbinitNcFile, Has_Header, Has_Structure, NotebookWriter):
xlims: Set the data limits for the x-axis in eV. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: Legend and title fontsize.
Returns:
matplotlib figure.
@ -182,12 +183,12 @@ class ScrFile(AbinitNcFile, Has_Header, Has_Structure, NotebookWriter):
set_axlims(ax, xlims, "x")
ax.grid(True)
ax.set_xlabel(r"$\omega$ [eV]")
ax.legend(loc="best")
ax.legend(loc="best", shadow=True, fontsize=fontsize)
return fig
@add_fig_kwargs
def plot_eelf(self, ax=None, xlims=None, **kwargs):
def plot_eelf(self, ax=None, xlims=None, fontsize=12, **kwargs):
r"""
Plot electron energy loss function.
@ -195,6 +196,7 @@ class ScrFile(AbinitNcFile, Has_Header, Has_Structure, NotebookWriter):
ax: matplotlib :class:`Axes` or None if a new figure should be created.
xlims: Set the data limits for the x-axis in eV. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used
fontsize: Legend and label fontsize:
Returns:
matplotlib figure.
@ -209,7 +211,7 @@ class ScrFile(AbinitNcFile, Has_Header, Has_Structure, NotebookWriter):
set_axlims(ax, xlims, "x")
ax.grid(True)
ax.set_xlabel(r"$\omega$ [eV]")
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@ -555,7 +557,7 @@ class _AwggMatrix(object):
return _latex_symbol_cplxmode(self.latex_name, cplx_mode)
@add_fig_kwargs
def plot_freq(self, gvec1, gvec2=None, waxis="real", cplx_mode="re-im", ax=None, **kwargs):
def plot_freq(self, gvec1, gvec2=None, waxis="real", cplx_mode="re-im", ax=None, fontsize=12, **kwargs):
r"""
Plot the frequency dependence of :math:`W_{G1, G2}(\omega)`
@ -568,6 +570,7 @@ class _AwggMatrix(object):
"angle" will display the phase of the complex number in radians.
Options can be concatenated with "-" e.g. "re-im"
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: legend and label fontsize.
Returns:
matplotlib figure.
@ -602,8 +605,8 @@ class _AwggMatrix(object):
ax.grid(True)
ax.set_xlabel(r"$\omega$ [eV]")
ax.set_title("%s, kpoint: %s" % (self.netcdf_name, self.kpoint))
ax.legend(loc="best")
ax.set_title("%s, kpoint: %s" % (self.netcdf_name, self.kpoint), fontsize=fontsize)
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig

1
abipy/electrons/tests/test_bse.py
View File

@ -118,6 +118,7 @@ class MdfRobotTest(AbipyTest):
plotter = robot.get_multimdf_plotter()
if self.has_matplotlib():
assert plotter.plot(show=False)
assert robot.plot(show=False)
#robot.plot_conv_mdf(self, hue, mdf_type="exc_mdf", **kwargs):
if self.has_nbformat():

5
abipy/electrons/tests/test_ebands.py
View File

@ -112,8 +112,12 @@ class ElectronBandsTest(AbipyTest):
self.assertMSONable(ni_ebands_kpath, test_if_subclass=False)
assert len(ni_ebands_kpath.to_json())
od = ni_ebands_kmesh.get_dict4pandas(with_spglib=False)
assert od["nsppol"] == 2 and od["nspinor"] == 1 and od["nspden"] == 2
df = ni_ebands_kpath.get_dataframe()
ni_ebands_kpath.to_xmgrace(self.get_tmpname(text=True))
ni_ebands_kpath.to_xmgrace(sys.stdout)
# BXSF cannot be produced because.
#ngkpt 6 6 6
@ -234,6 +238,7 @@ class ElectronBandsTest(AbipyTest):
si_edos = si_ebands_kmesh.get_edos()
repr(si_edos); str(si_edos)
assert ElectronDos.as_edos(si_edos, {}) is si_edos
assert si_edos == si_edos and not (si_edos != si_edos)
edos_samevals = ElectronDos.as_edos(si_ebands_kmesh, {})
assert ElectronDos.as_edos(si_ebands_kmesh, {}) == si_edos
assert ElectronDos.as_edos(abidata.ref_file("si_scf_GSR.nc"), {}) == si_edos

38
abipy/electrons/tests/test_gsr.py
View File

@ -142,7 +142,7 @@ class GstRobotTest(AbipyTest):
gsr_path = abidata.ref_file("si_scf_GSR.nc")
robot = abilab.GsrRobot()
robot.add_file("gsr0", gsr_path)
assert len(robot.ncfiles) == 1
assert len(robot.abifiles) == 1
assert robot.EXT == "GSR"
repr(robot); str(robot)
@ -169,7 +169,15 @@ class GstRobotTest(AbipyTest):
if self.has_matplotlib():
assert ebands_plotter.gridplot(show=False)
assert robot.combiplot_ebands(show=False)
assert robot.gridplot_ebands(show=False)
assert robot.boxplot_ebands(show=False)
assert robot.combiboxplot_ebands(show=False)
assert edos_plotter.gridplot(show=False)
assert robot.combiplot_edos(show=False)
assert robot.gridplot_edos(show=False)
# Get pandas dataframe.
df = robot.get_dataframe()
@ -177,9 +185,33 @@ class GstRobotTest(AbipyTest):
self.assert_equal(df["ecut"].values, 6.0)
self.assert_almost_equal(df["energy"].values, -241.2364683)
# FIXME
#eos = robot.eos_fit()
if self.has_matplotlib():
assert robot.plot_xy_with_hue(df, x="nkpt", y="pressure", hue="a", show=False)
# Note: This is not a real EOS since we have a single volume.
# But testing is better than not testing.
r = robot.get_eos_fits_dataframe()
assert hasattr(r, "fits") and hasattr(r, "dataframe")
if self.has_matplotlib():
assert robot.gridplot_eos(show=False)
if self.has_nbformat():
robot.write_notebook(nbpath=self.get_tmpname(text=True))
robot.close()
# Test other class methods
filepath = abidata.ref_file("si_scf_GSR.nc")
robot = abilab.GsrRobot.from_dirs(os.path.dirname(filepath), abspath=True)
assert len(robot) == 2
assert robot[filepath].filepath == filepath
# Test from_dir_glob
pattern = "%s/*/si_ebands/" % abidata.dirpath
same_robot = abilab.GsrRobot.from_dir_glob(pattern, abspath=True)
assert len(same_robot) == 2
assert set(robot.labels) == set(same_robot.labels)
robot.close()
same_robot.close()

7
abipy/electrons/tests/test_gw.py
View File

@ -81,6 +81,7 @@ class TestSigresFile(AbipyTest):
"""Test SIGRES File."""
sigres = abilab.abiopen(abidata.ref_file("tgw1_9o_DS4_SIGRES.nc"))
assert sigres.nsppol == 1
sigres.print_qps(precision=5, ignore_imag=False)
# In this run IBZ = kptgw
assert len(sigres.ibz) == 6
@ -99,8 +100,8 @@ class TestSigresFile(AbipyTest):
qpgaps = [3.53719151871085, 4.35685250045637, 4.11717896881632,
8.71122659251508, 3.29693118466282, 3.125545059031]
self.assert_almost_equal(sigres.qpgaps, np.reshape(qpgaps, (1, 6)))
ik = 2
df = sigres.get_dataframe_sk(spin=0, kpoint=ik)
same_df = sigres.get_dataframe_sk(spin=0, kpoint=sigres.gwkpoints[ik])
@ -236,10 +237,6 @@ class SigresRobotTest(AbipyTest):
df_sk = robot.merge_dataframes_sk(spin=0, kpoint=[0, 0, 0])
qpdata = robot.get_qpgaps_dataframe(with_geo=True)
# This cannot be tested because it changes the matplotlib backend
if self.has_seaborn():
robot.plot_conv_qpgap(x_vars="sigma_nband", show=False)
if self.has_nbformat():
robot.write_notebook(nbpath=self.get_tmpname(text=True))

12
abipy/eph/eph.py
View File

@ -183,7 +183,7 @@ class A2f(object):
@add_fig_kwargs
def plot(self, units="eV", with_lambda=True, exchange_xy=False, ax=None,
xlims=None, ylims=None, label=None, **kwargs):
xlims=None, ylims=None, label=None, fontsize=12, **kwargs):
"""
Plot a2F(w), it's primitive lambda(w) and optionally the individual
contributions due to the phonon branches.
@ -198,6 +198,7 @@ class A2f(object):
or scalar e.g. `left`. If left (right) is None, default values are used
ylims: Limits for y-axis. See xlims for API.
label: True to add legend label to each curve.
fontsize: Legend and title fontsize
Returns:
`matplotlib` figure
@ -243,7 +244,7 @@ class A2f(object):
ax.grid(True)
set_axlims(ax, xlims, "x")
set_axlims(ax, ylims, "y")
if label: ax.legend(loc="best")
if label: ax.legend(loc="best", shadow=True, fontsize=fontsize)
return fig
@ -561,7 +562,7 @@ class EphFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter):
self.reader.close()
@add_fig_kwargs
def plot_eph_strength(self, what="lambda", ylims=None, ax=None, label=None, **kwargs):
def plot_eph_strength(self, what="lambda", ylims=None, ax=None, label=None, fontsize=12, **kwargs):
"""
Plot phonon bands with eph coupling strength lambda(q, nu)
@ -571,6 +572,7 @@ class EphFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter):
or scalar e.g. `left`. If left (right) is None, default values are used
ax: matplotlib :class:`Axes` or None if a new figure should be created.
label: String used to label the plot in the legend.
fontsize: Legend and title fontsize.
Returns:
`matplotlib` figure
@ -607,7 +609,7 @@ class EphFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter):
ax.set_ylabel(ylabel)
set_axlims(ax, ylims, "y")
if label: ax.legend(loc="best")
if label: ax.legend(loc="best", shadow=True, fontsize=fontsize)
return fig
@ -958,7 +960,7 @@ class EphRobot(Robot, RobotWithEbands, RobotWithPhbands):
nbv.new_code_cell("robot.plot_a2f_convergence();"),
])
if all(ncf.has_a2ftr for ncf in self.ncfiles):
if all(ncf.has_a2ftr for ncf in self.abifiles):
nb.cells.extend([
nbv.new_code_cell("robot.plot_a2ftr_convergence();"),
])

31
abipy/eph/sigeph.py
View File

@ -13,7 +13,6 @@ import pandas as pd
import pymatgen.core.units as units
import abipy.core.abinit_units as abu
from collections import OrderedDict, namedtuple
from monty.string import marquee, list_strings
from monty.functools import lazy_property
@ -129,7 +128,7 @@ class QpTempState(namedtuple("QpTempState", "tmesh e0 qpe ze0 spin kpoint band")
num_plots, ncols, nrows = len(fields), 1, 1
if num_plots > 1:
ncols = 2
nrows = (num_plots//ncols) + (num_plots % ncols)
nrows = (num_plots // ncols) + (num_plots % ncols)
# Build grid of plots.
ax_list, fig, plt = get_axarray_fig_plt(ax_list, nrows=nrows, ncols=ncols,
@ -158,7 +157,7 @@ class QpTempState(namedtuple("QpTempState", "tmesh e0 qpe ze0 spin kpoint band")
# Get around a bug in matplotlib
if num_plots % ncols != 0:
ax_list[-1].plot([0,1], [0,1], lw=0)
ax_list[-1].plot([0, 1], [0, 1], lw=0)
ax_list[-1].axis('off')
if label is not None:
@ -304,7 +303,7 @@ class QpTempList(list):
num_plots, ncols, nrows = len(fields), 1, 1
if num_plots > 1:
ncols = 2
nrows = (num_plots//ncols) + (num_plots % ncols)
nrows = (num_plots // ncols) + (num_plots % ncols)
# Build grid of plots.
ax_list, fig, plt = get_axarray_fig_plt(ax_list, nrows=nrows, ncols=ncols, sharex=True, sharey=False, squeeze=False)
@ -336,7 +335,7 @@ class QpTempList(list):
# Get around a bug in matplotlib
if num_plots % ncols != 0:
ax_list[-1].plot([0,1], [0,1], lw=0)
ax_list[-1].plot([0, 1], [0, 1], lw=0)
ax_list[-1].axis('off')
#if label is not None:
@ -448,7 +447,7 @@ class EphSelfEnergy(object):
@add_fig_kwargs
def plot_tdep(self, itemps="all", zero_energy="e0", cmap="jet", ax_list=None,
what_list=("re", "im", "spfunc"), xlims=None, **kwargs):
what_list=("re", "im", "spfunc"), xlims=None, fontsize=12, **kwargs):
"""
Plot the real/imaginary part of self-energy as well as the spectral function for
the different temperatures with a color map.
@ -461,6 +460,7 @@ class EphSelfEnergy(object):
what_list:
xlims: Set the data limits for the x-axis. Accept tuple e.g. `(left, right)`
or scalar e.g. `left`. If left (right) is None, default values are used.
fontsize: legend and label fontsize.
Returns:
`matplotlib` figure
@ -480,7 +480,7 @@ class EphSelfEnergy(object):
color=cmap(itemp / self.ntemp),
label=tlabels[itemp] if i == 0 else None,
)
if i == 0: ax.legend(loc="best")
if i == 0: ax.legend(loc="best", shadow=True, fontsize=fontsize)
set_axlims(ax, xlims, "x")
return fig
@ -665,12 +665,13 @@ class SigEPhFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter)
#def get_dirgaps_dataframe(self):
@add_fig_kwargs
def plot_qpgaps_t(self, ax=None, **kwargs):
def plot_qpgaps_t(self, ax=None, fontsize=12, **kwargs):
"""
Plot the KS and the QP(T) direct gaps for all the k-points available on file.
Args:
ax: matplotlib :class:`Axes` or None if a new figure should be created.
fontsize: legend and title fontsize.
Returns:
`matplotlib` figure
@ -695,7 +696,7 @@ class SigEPhFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter)
ax.grid(True)
ax.set_xlabel("Temperature [K]")
ax.set_ylabel("Direct gap [eV]")
ax.legend(loc="best")
ax.legend(loc="best", fontsize=fontsize, shadow=True)
return fig
@ -756,16 +757,16 @@ class SigEPhRobot(Robot, RobotWithEbands):
def __init__(self, *args):
super(SigEPhRobot, self).__init__(*args)
if len(self.ncfiles) in (0, 1): return
if len(self.abifiles) in (0, 1): return
# Check dimensions and self-energy states and issue warning.
warns = []; wapp = warns.append
nc0 = self.ncfiles[0]
nc0 = self.abifiles[0]
same_nsppol, same_nkcalc = True, True
if any(nc.nsppol != nc0.nsppol for nc in self.ncfiles):
if any(nc.nsppol != nc0.nsppol for nc in self.abifiles):
same_nsppol = False
wapp("Comparing ncfiles with different values of nsppol.")
if any(nc.nkcalc != nc0.nkcalc for nc in self.ncfiles):
if any(nc.nkcalc != nc0.nkcalc for nc in self.abifiles):
same_nkcalc = False
wapp("Comparing ncfiles with different number of k-points in self-energy. Do")
@ -774,9 +775,9 @@ class SigEPhRobot(Robot, RobotWithEbands):
# Different values of bstart_ks are difficult to handle
# Because the high-level API assumes an absolute global index
# Should decided how to treat this case: either raise or interpret band as an absolute band index.
if any(np.any(nc.bstart_sk != nc0.bstart_sk) for nc in self.ncfiles):
if any(np.any(nc.bstart_sk != nc0.bstart_sk) for nc in self.abifiles):
wapp("Comparing ncfiles with different values of bstart_sk")
if any(np.any(nc.bstop_sk != nc0.bstop_sk) for nc in self.ncfiles):
if any(np.any(nc.bstop_sk != nc0.bstop_sk) for nc in self.abifiles):
wapp("Comparing ncfiles with different values of bstop_sk")
if wapp:

15
abipy/examples/flows/README.txt
View File

@ -4,6 +4,15 @@ AbiPy Flow Gallery
==================
This gallery contains python scripts to generate AbiPy flows from the command line.
Run the script to generate the directory with the flow.
Use `-s` to generate the flow and run it with the scheduler.
Use `--help` for see available options.
Run the scripts to generate the directory with the flow and then
use ``abirun.py`` to execute the flow.
Alternatively, one can use the `-s` option to generate the flow and run it immediately with the scheduler.
Use `--help` for available options.
WARNING:
The following examples show how to use python and the AbiPy API to generate and run
Abinit calculations in a semi-automatic way.
These examples are not supposed to produce physically meaningful results
as input parameters are usually underconverged.

46
abipy/examples/flows/develop/run_eph_al.py
View File

@ -27,7 +27,7 @@ def build_flow(options):
pseudos = abidata.pseudos("Al.oncvpsp")
structure = abilab.Structure.from_abivars(
acell=3*[7.5],
acell=3 * [7.5],
rprim=[0.0, 0.5, 0.5,
0.5, 0.0, 0.5,
0.5, 0.5, 0.0],
@ -37,6 +37,9 @@ def build_flow(options):
znucl=13,
)
same_structure = abilab.Structure.fcc(a=7.5, species=["Al"], units="bohr")
assert same_structure == structure
gs_inp = abilab.AbinitInput(structure, pseudos)
gs_inp.set_vars(
@ -66,31 +69,30 @@ def build_flow(options):
ph_work = flowtk.PhononWork.from_scf_task(work0[0], qpoints)
flow.register_work(ph_work)
# Build input file for E-PH run.
eph_inp = gs_inp.new_with_vars(
optdriver=7,
ddb_ngqpt=ddb_ngqpt, # q-mesh used to produce the DDB file (must be consistent with DDB data)
eph_intmeth=2, # Tetra method
eph_fsewin="0.8 eV", # Energy window around Ef
eph_mustar=0.12, # mustar parameter
)
# Set q-path for phonons and phonon linewidths.
eph_inp.set_qpath(20)
# TODO: Define wstep and smear
# Set q-mesh for phonons DOS and a2F(w)
eph_inp.set_phdos_qmesh(nqsmall=16, method="tetra")
eph_work = flow.register_work(flowtk.Work())
eph_deps = {work0[0]: "WFK", ph_work: ["DDB", "DVDB"]}
eph_task = eph_work.register_eph_task(eph_inp, deps=eph_deps)
flow.allocate()
flow.use_smartio()
for eph_ngqpt_fine in [(4, 4, 4,), (8, 8, 8)]:
# Build input file for E-PH run.
eph_inp = gs_inp.new_with_vars(
optdriver=7,
ddb_ngqpt=ddb_ngqpt, # q-mesh used to produce the DDB file (must be consistent with DDB data)
eph_intmeth=2, # Tetra method
eph_fsewin="0.8 eV", # Energy window around Ef
eph_mustar=0.12, # mustar parameter
eph_ngqpt_fine=eph_ngqpt_fine, # Interpolate DFPT potentials if != ddb_ngqpt
)
# Set q-path for phonons and phonon linewidths.
eph_inp.set_qpath(20)
# TODO: Define wstep and smear
# Set q-mesh for phonons DOS and a2F(w)
eph_inp.set_phdos_qmesh(nqsmall=16, method="tetra")
eph_work.register_eph_task(eph_inp, deps=eph_deps)
flow.allocate(use_smartio=True)
# EPH does not support autoparal (yet)
#eph_task.with_fixed_mpi_omp(1, 1)
return flow

5
abipy/examples/flows/develop/run_gruneisen_phonopy_si.py
View File

@ -31,14 +31,15 @@ def build_flow(options):
# Build input for GS calculation.
gsinp = abilab.AbinitInput(structure, pseudos)
gsinp.set_vars(ecut=8, nband=4, toldff=1.e-6)
# This gives ngkpt = 4x4x4 with 4 shifts for the initial unit cell.
# The k-point sampling will be rescaled when we build the supercell in PhonopyWork.
gsinp.set_autokmesh(nksmall=4)
#gsinp.set_vars(ngkpt=[4, 4, 4])
flow = flowtk.Flow(workdir=options.workdir)
# Use a 2x2x2 supercell to compute phonons with phonopy
work = PhonopyGruneisenWork.from_gs_input(gsinp, voldelta=0.01, scdims=[4, 4, 4])
work = PhonopyGruneisenWork.from_gs_input(gsinp, voldelta=0.01, scdims=[2, 2, 2])
flow.register_work(work)
return flow

151
abipy/examples/flows/develop/run_optic_old.py Executable file
View File

@ -0,0 +1,151 @@
#!/usr/bin/env python
r"""
Optic Flow
==========
Optical spectra with Optic.
"""
from __future__ import print_function, division, unicode_literals, absolute_import
import sys
import os
import abipy.data as abidata
import abipy.abilab as abilab
import abipy.flowtk as flowtk
def build_flow(options, paral_kgb=0):
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
multi = abilab.MultiDataset(structure=abidata.structure_from_ucell("GaAs"),
pseudos=abidata.pseudos("31ga.pspnc", "33as.pspnc"), ndtset=5)
# Global variables
kmesh = dict(ngkpt=[4, 4, 4],
nshiftk=4,
shiftk=[[0.5, 0.5, 0.5],
[0.5, 0.0, 0.0],
[0.0, 0.5, 0.0],
[0.0, 0.0, 0.5]]
)
global_vars = dict(ecut=2, paral_kgb=paral_kgb)
global_vars.update(kmesh)
multi.set_vars(global_vars)
# Dataset 1 (GS run)
multi[0].set_vars(
tolvrs=1e-6,
nband=4,
)
# NSCF run with large number of bands, and points in the the full BZ
multi[1].set_vars(
iscf=-2,
nband=20,
nstep=25,
kptopt=1,
tolwfr=1.e-9,
#kptopt=3,
)
# Fourth dataset : ddk response function along axis 1
# Fifth dataset : ddk response function along axis 2
# Sixth dataset : ddk response function along axis 3
for idir in range(3):
rfdir = 3 * [0]
rfdir[idir] = 1
multi[2+idir].set_vars(
iscf=-3,
nband=20,
nstep=1,
nline=0,
prtwf=3,
kptopt=3,
nqpt=1,
qpt=[0.0, 0.0, 0.0],
rfdir=rfdir,
rfelfd=2,
tolwfr=1.e-9,
)
scf_inp, nscf_inp, ddk1, ddk2, ddk3 = multi.split_datasets()
# Initialize the flow.
flow = flowtk.Flow(options.workdir, manager=options.manager)
bands_work = flowtk.BandStructureWork(scf_inp, nscf_inp)
flow.register_work(bands_work)
ddk_work = flowtk.Work()
for inp in [ddk1, ddk2, ddk3]:
ddk_work.register_ddk_task(inp, deps={bands_work.nscf_task: "WFK"})
flow.register_work(ddk_work)
# Optic does not support MPI with ncpus > 1.
optic_input = abilab.OpticInput(
broadening=0.002, # Value of the smearing factor, in Hartree
domega=0.0003, # Frequency mesh.
maxomega=0.3,
scissor=0.000, # Scissor shift if needed, in Hartree
tolerance=0.002, # Tolerance on closeness of singularities (in Hartree)
num_lin_comp=1, # Number of components of linear optic tensor to be computed
lin_comp=11, # Linear coefficients to be computed (x=1, y=2, z=3)
num_nonlin_comp=2, # Number of components of nonlinear optic tensor to be computed
nonlin_comp=(123, 222), # Non-linear coefficients to be computed
)
# TODO
# Check is the order of the 1WF files is relevant. Can we use DDK files ordered
# in an arbitrary way or do we have to pass (x,y,z)?
optic_task = flowtk.OpticTask(optic_input, nscf_node=bands_work.nscf_task, ddk_nodes=ddk_work)
flow.register_task(optic_task)
return flow
def optic_flow_from_files():
# Optic does not support MPI with ncpus > 1.
manager = flowtk.TaskManager.from_user_config()
manager.set_mpi_procs(1)
flow = flowtk.Flow(workdir="OPTIC_FROM_FILE", manager=manager)
ddk_nodes = [
"/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_1/task_0/outdata/out_1WF",
"/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_1/task_1/outdata/out_1WF",
"/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_1/task_2/outdata/out_1WF",
]
nscf_node = "/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_0/task_1/outdata/out_WFK"
optic_task = flowtk.OpticTask(optic_input, nscf_node=nscf_node, ddk_nodes=ddk_nodes)
flow.register_task(optic_task)
return flow
# This block generates the thumbnails in the Abipy gallery.
# You can safely REMOVE this part if you are using this script for production runs.
if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
@flowtk.flow_main
def main(options):
"""
This is our main function that will be invoked by the script.
flow_main is a decorator implementing the command line interface.
Command line args are stored in `options`.
"""
return build_flow(options)
if __name__ == "__main__":
sys.exit(main())

2
abipy/examples/flows/run_qptdmscr.py → abipy/examples/flows/develop/run_qptdmscr.py
View File

@ -130,7 +130,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
#import tempfile
#options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
#build_flow(options).plot_networkx(tight_layout=True)
#build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main

36
abipy/examples/flows/run_wyckoff.py → abipy/examples/flows/develop/run_wyckoff.py
View File

@ -4,7 +4,7 @@ Wyckoff Flow
============
This example shows how to compute the band structure of a set of
crystalline structures obtained by changing a set of internal paramaters
crystalline structures obtained by changing the set of internal parameters (wyckoff positions).
"""
from __future__ import division, print_function, unicode_literals, absolute_import
@ -19,14 +19,13 @@ def special_positions(lattice, u):
"""Construct the crystalline `Structure` for given value of the internal parameter u."""
frac_coords = {}
frac_coords["Si"] = [ [1/2 + u, 1/2 - u, 0],
[u, -u, u],
]
frac_coords["Si"] = [
[1/2 + u, 1/2 - u, 0],
[u, -u, u],
]
frac_coords["O"] = [ [0, 0, u] ]
species, coords = [], []
for symbol, positions in frac_coords.items():
species += len(positions) * [symbol]
coords += positions
@ -37,13 +36,12 @@ def special_positions(lattice, u):
def build_flow(options):
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_","flow_")
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
pseudos = abidata.pseudos("14si.pspnc", "8o.pspnc")
base_structure = abilab.Structure.from_file(abidata.cif_file("si.cif"))
news, uparams = [], [0.2, 0.3]
for u in uparams:
new = special_positions(base_structure.lattice, u)
news.append(new)
@ -63,32 +61,30 @@ def make_workflow(structure, pseudos, paral_kgb=1):
Return a `Workflow` object defining a band structure calculation
for given `Structure`.
"""
# Variables global to the SCF and the NSCF run.
global_vars = dict(
ecut=abilab.FloatWithUnit(100, "eV").to("Ha"),
paral_kgb=paral_kgb,
#nband=8,
)
# GS + NSCF run
multi = abilab.MultiDataset(structure, pseudos=pseudos, ndtset=2)
multi.set_vars(global_vars)
nval = structure.num_valence_electrons(pseudos)
# Variables global to the SCF and the NSCF run.
multi.set_vars(
ecut=15,
paral_kgb=paral_kgb,
nband=nval//2 + 4, # occupied + 4 empty
)
# (GS run)
multi[0].set_kmesh(ngkpt=[8,8,8], shiftk=[0,0,0])
multi[0].set_kmesh(ngkpt=[4, 4, 4], shiftk=[0, 0, 0])
multi[0].set_vars(tolvrs=1e-6)
# (NSCF run)
multi[1].set_vars(
iscf=-2,
tolwfr=1e-12,
kptopt=0,
nkpt=1,
kpt=[0, 0, 0],
)
multi[1].set_kpath(ndivsm=8)
gs_inp, nscf_inp = multi.split_datasets()
return flowtk.BandStructureWork(gs_inp, nscf_inp)

155
abipy/examples/flows/run_eos.py Executable file
View File

@ -0,0 +1,155 @@
#!/usr/bin/env python
r"""
Flow for Equation of State
==========================
Flow to compute the equation of state by fitting E(V) at T = 0.
"""
from __future__ import division, print_function, unicode_literals, absolute_import
import sys
import os
import abipy.data as abidata
import abipy.abilab as abilab
import abipy.flowtk as flowtk
def build_flow(options):
# Set working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
# Build GS input file.
pseudos = abidata.pseudos("Si.GGA_PBE-JTH-paw.xml")
#silicon = abilab.Structure.zincblende(5.431, ["Si", "Si"], units="ang")
silicon = abidata.cif_file("si.cif")
scf_input = abilab.AbinitInput(silicon, pseudos)
ecut = 12
scf_input.set_vars(
ecut=ecut,
pawecutdg=40,
nband=6,
paral_kgb=0,
iomode=3,
toldfe=1e-9,
)
# K-point sampling (shifted)
scf_input.set_autokmesh(nksmall=4)
from abipy.flowtk.gs_works import EosWork
flow = flowtk.Flow(options.workdir, manager=options.manager)
# Si is cubic and atomic positions are fixed by symmetry so we
# use move_atoms=False to compute E(V) with SCF-GS tasks instead of
# performing a constant-volume optimization of the cell geometry.
work = EosWork.from_scf_input(scf_input, move_atoms=False, ecutsm=0.5)
flow.register_work(work)
flow.allocate(use_smartio=True)
return flow
# This block generates the thumbnails in the Abipy gallery.
# You can safely REMOVE this part if you are using this script for production runs.
if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
def main(options):
"""
This is our main function that will be invoked by the script.
flow_main is a decorator implementing the command line interface.
Command line args are stored in `options`.
"""
return build_flow(options)
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_eos.py -s
#
# then use:
#
# abirun.py flow_eos structures
#
# to get info about the intial and final structures.
#
# .. code-block:: bash
#
# Lattice parameters:
# formula natom angle0 angle1 angle2 a b c volume \
# w0_t0_in Si2 2 60.0 60.0 60.0 3.854 3.854 3.854 40.479
# w0_t0_out Si2 2 60.0 60.0 60.0 3.854 3.854 3.854 40.479
# w0_t1_in Si2 2 60.0 60.0 60.0 3.857 3.857 3.857 40.582
# w0_t1_out Si2 2 60.0 60.0 60.0 3.857 3.857 3.857 40.582
# w0_t2_in Si2 2 60.0 60.0 60.0 3.861 3.861 3.861 40.684
# w0_t2_out Si2 2 60.0 60.0 60.0 3.861 3.861 3.861 40.684
# w0_t3_in Si2 2 60.0 60.0 60.0 3.864 3.864 3.864 40.786
# w0_t3_out Si2 2 60.0 60.0 60.0 3.864 3.864 3.864 40.786
# w0_t4_in Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w0_t4_out Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w0_t5_in Si2 2 60.0 60.0 60.0 3.870 3.870 3.870 40.991
# w0_t5_out Si2 2 60.0 60.0 60.0 3.870 3.870 3.870 40.991
# w0_t6_in Si2 2 60.0 60.0 60.0 3.873 3.873 3.873 41.093
# w0_t6_out Si2 2 60.0 60.0 60.0 3.873 3.873 3.873 41.093
# w0_t7_in Si2 2 60.0 60.0 60.0 3.877 3.877 3.877 41.195
# w0_t7_out Si2 2 60.0 60.0 60.0 3.877 3.877 3.877 41.195
# w0_t8_in Si2 2 60.0 60.0 60.0 3.880 3.880 3.880 41.297
# w0_t8_out Si2 2 60.0 60.0 60.0 3.880 3.880 3.880 41.297
#
# abispg_num P [GPa] Max|F| eV/ang task_class status
# w0_t0_in None NaN NaN ScfTask Completed
# w0_t0_out 227 1.327 9.506e-28 ScfTask Completed
# w0_t1_in None NaN NaN ScfTask Completed
# w0_t1_out 227 1.092 1.877e-27 ScfTask Completed
# w0_t2_in None NaN NaN ScfTask Completed
# w0_t2_out 227 0.859 4.649e-27 ScfTask Completed
# w0_t3_in None NaN NaN ScfTask Completed
# w0_t3_out 227 0.629 5.062e-27 ScfTask Completed
# w0_t4_in None NaN NaN ScfTask Completed
# w0_t4_out 227 0.403 5.945e-27 ScfTask Completed
# w0_t5_in None NaN NaN ScfTask Completed
# w0_t5_out 227 0.178 6.306e-27 ScfTask Completed
# w0_t6_in None NaN NaN ScfTask Completed
# w0_t6_out 227 -0.042 1.738e-27 ScfTask Completed
# w0_t7_in None NaN NaN ScfTask Completed
# w0_t7_out 227 -0.259 2.725e-27 ScfTask Completed
# w0_t8_in None NaN NaN ScfTask Completed
# w0_t8_out 227 -0.474 6.939e-27 ScfTask Completed
#
# Use `--verbose` to print atoms.
#
# Silicon is a cubic materials and the positions are fixed by symmeetry so we only need ScfTask to get E(V).
#
# The EOS results are saved in the JSON file in flow_eos/w0/outdata/eos_data.json.
#
# We can also use the GSR robot to compute and plot the EOS.
#
# Use
#
# abirun.py flow_eos robot GSR
#
# to build a robot for GSR files and start an ipython shell.
# Then, inside ipython, type:
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: robot.gridplot_eos()
#
# to produce
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_eos.png?raw=true
# :alt: Equation of state of Si obtained with different models.
#

18
abipy/examples/flows/run_fe_ebands.py
View File

@ -74,7 +74,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -89,3 +89,19 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_fe_ebands -s
#
# then use:
#
# abirun.py flow_fe_ebands ebands -p
#
# to analyze all the band structures produced by the Flow and plot the data
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_fe_ebands.png?raw=true
# :alt: Band structures of Fe with nsppol 1, 2.
#

134
abipy/examples/flows/run_gaas_ebands_soc.py Executable file
View File

@ -0,0 +1,134 @@
#!/usr/bin/env python
r"""
Band structure with/without spin-orbit
======================================
This example shows how to compute the band structure of GaAs with and without spin-orbit term.
We essentially build two BandStructureWork inside a loop over nspinor in [1, 2]
nspinor = 1 corresponds to a standard collinear calculation for non-magnetic systems while
nspinor = 2 gives us the non-collinear case with spinor wavefunctions required for the treatment of SOC.
Some of the variables in the input files must be changed depending on the value of nspinor.
We use relativistic NC pseudos made of two terms: scalar pseudo + SOC term.
The SOC term can be deactivated with the input variable `so_psp`.
"""
from __future__ import division, print_function, unicode_literals, absolute_import
import sys
import os
import abipy.data as abidata
import abipy.abilab as abilab
import abipy.flowtk as flowtk
def build_flow(options):
# Set working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
structure = abidata.structure_from_ucell("GaAs")
pseudos = abidata.pseudos("Ga-low_r.psp8", "As_r.psp8")
num_electrons = structure.num_valence_electrons(pseudos)
#print("num_electrons:", num_electrons)
# Usa same shifts in all tasks.
ngkpt = [4, 4, 4]
shiftk= [
[0.5, 0.5, 0.5],
[0.5, 0.0, 0.0],
[0.0, 0.5, 0.0],
[0.0, 0.0, 0.5],
]
# NSCF run on k-path with large number of bands
kptbounds = [
[0.5, 0.0, 0.0], # L point
[0.0, 0.0, 0.0], # Gamma point
[0.0, 0.5, 0.5], # X point
]
# Initialize the flow.
flow = flowtk.Flow(options.workdir, manager=options.manager)
for nspinor in [1, 2]:
# Multi will contain two datasets (GS + NSCF) for the given nspinor.
multi = abilab.MultiDataset(structure=structure, pseudos=pseudos, ndtset=2)
# Global variables.
multi.set_vars(
ecut=20,
nspinor=nspinor,
nspden=1 if nspinor == 1 else 4,
so_psp="*0" if nspinor == 1 else "*1", # Important!
#paral_kgb=1,
)
nband_occ = num_electrons // 2 if nspinor == 1 else num_electrons
#print(nband_occ)
# Dataset 1 (GS run)
multi[0].set_vars(tolvrs=1e-8, nband=nband_occ + 4)
multi[0].set_kmesh(ngkpt=ngkpt, shiftk=shiftk, kptopt=1 if nspinor == 1 else 4)
multi[1].set_vars(iscf=-2, nband=nband_occ + 4, tolwfr=1.e-12)
multi[1].set_kpath(ndivsm=10, kptbounds=kptbounds)
# Get the SCF and the NSCF input.
scf_input, nscf_input = multi.split_datasets()
flow.register_work(flowtk.BandStructureWork(scf_input, nscf_input))
return flow
# This block generates the thumbnails in the Abipy gallery.
# You can safely REMOVE this part if you are using this script for production runs.
if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
def main(options):
"""
This is our main function that will be invoked by the script.
flow_main is a decorator implementing the command line interface.
Command line args are stored in `options`.
"""
return build_flow(options)
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_gaas_ebands_soc.py -s
#
# then use:
#
# abirun.py flow_si_ebands ebands --plot -t NscfTask
#
# to analyze (and plot) the electronic bands produced by the NsfTasks of the Flow.
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_gaas_ebands_soc.png?raw=true
# :alt: Band structure of GaAs without/with SOC.
#
# Alternatively, one can start a GSR robot for the NscfTask with:
#
# abirun.py flow_gaas_ebands_soc/ robot GSR -t NscfTask
#
# and then plot the two band structure on the same figure with:
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: robot.combiplot_ebands()
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_gaas_ebands_soc_combiplot.png?raw=true
# :alt: Band structure of GaAs without/with SOC on the same graph

6
abipy/examples/flows/run_gwconv_ecuteps.py
View File

@ -83,7 +83,7 @@ def make_inputs(paral_kgb=1):
def build_flow(options):
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_","flow_")
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
# Get our templates
scf_inp, nscf_inp, scr_inp, sig_inp = make_inputs()
@ -120,14 +120,13 @@ def build_flow(options):
return flow
# This block generates the thumbnails in the Abipy gallery.
# You can safely REMOVE this part if you are using this script for production runs.
if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -139,6 +138,5 @@ def main(options):
"""
return build_flow(options)
if __name__=="__main__":
sys.exit(main())

72
abipy/examples/flows/run_ht_si_bsemdf.py
View File

@ -3,7 +3,7 @@ r"""
Bethe-Salpeter Flow with factory functions
==========================================
Calculation of the BSE spectrum with the HT interface.
Calculation of the BSE spectrum with the High-throuhput interface (factory functions).
"""
from __future__ import print_function, division, unicode_literals, absolute_import
@ -26,11 +26,11 @@ def build_flow(options):
kppa = scf_kppa = 1
nscf_nband = 6
nscf_ngkpt = [4,4,4]
nscf_ngkpt = [4, 4, 4]
nscf_shiftk = [0.1, 0.2, 0.3]
bs_loband = 2
bs_nband = nscf_nband
mbpt_sciss = 0.7
mbpt_sciss = 0.7 * abilab.units.eV_to_Ha
mdf_epsinf = 12
ecuteps = 2
ecut = 12
@ -42,19 +42,13 @@ def build_flow(options):
structure, pseudos,
scf_kppa, nscf_nband, nscf_ngkpt, nscf_shiftk,
ecuteps, bs_loband, bs_nband, mbpt_sciss, mdf_epsinf,
ecut=ecut,#$ pawecutdg=None,
ecut=ecut,# pawecutdg=None,
exc_type="TDA", bs_algo="haydock", accuracy="normal", spin_mode="unpolarized",
smearing=None)
#smearing="fermi_dirac:0.1 eV", charge=0.0, scf_algorithm=None)
work = flowtk.BseMdfWork(scf_input=multi[0], nscf_input=multi[1], bse_inputs=multi[2:])
#from pymatgen.io.abinit.calculations import bse_with_mdf_work
#work = bse_with_mdf_work(structure, pseudos, scf_kppa, nscf_nband, nscf_ngkpt, nscf_shiftk,
# ecuteps, bs_loband, bs_nband, mbpt_sciss, mdf_epsinf,
# accuracy="normal", spin_mode="unpolarized", smearing=None,
# charge=0.0, scf_solver=None, **extra_abivars)
flow.register_work(work)
return flow
@ -65,8 +59,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -81,3 +74,58 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_ht_si_bsemdf.py -s
#
# then use:
#
# abirun.py flow_ht_si_bsemdf/ status -v
#
# to get the status of the flow
#
# .. code-block:: bash
#
# Work #0: <BseMdfWork, node_id=241906, workdir=flow_ht_si_bsemdf/w0>, Finalized=True
# +--------+-----------+-----------------+--------------+------------+----------+-----------------+----------+-----------+
# | Task | Status | Queue | MPI|Omp|Gb | Warn|Com | Class | Sub|Rest|Corr | Time | Node_ID |
# +========+===========+=================+==============+============+==========+=================+==========+===========+
# | w0_t0 | Completed | 62430@localhost | 1| 1|2.0 | 2| 0 | ScfTask | (1, 0, 0) | 0:00:01R | 241907 |
# +--------+-----------+-----------------+--------------+------------+----------+-----------------+----------+-----------+
# | w0_t1 | Completed | 62438@localhost | 2| 1|2.0 | 3| 0 | NscfTask | (1, 0, 0) | 0:00:05R | 241908 |
# +--------+-----------+-----------------+--------------+------------+----------+-----------------+----------+-----------+
# | w0_t2 | Completed | 62449@localhost | 2| 1|2.0 | 10| 6 | BseTask | (1, 0, 0) | 0:00:17R | 241909 |
# +--------+-----------+-----------------+--------------+------------+----------+-----------------+----------+-----------+
#
# all_ok reached
#
# The macroscopic dielectric function produced by the BseTask is stored in the ``out_MDF.nc`` file:
#
# .. code-block:: bash
#
# abirun.py flow_ht_si_bsemdf/ listext MDF
#
# Found 1 files with extension `MDF` produced by the flow
# File Size [Mb] Node_ID Node Class
# ------------------------------------------ ----------- --------- ------------
# flow_ht_si_bsemdf/w0/t2/outdata/out_MDF.nc 0.4 241909 BseTask
#
# Open the file with:
#
# abiopen.py flow_ht_si_bsemdf/w0/t2/outdata/out_MDF.nc
#
# and then inside the ipython terminal, issue:
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: abifile.plot_mdfs()
#
# to produce:
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_ht_si_bsemdf.png?raw=true
# :alt: Band structure of Si in the IBZ and along a k-path
#

55
abipy/examples/flows/run_ht_si_ebands.py
View File

@ -3,7 +3,8 @@ r"""
Band structure Flow with factory functions
==========================================
Band structure of silicon with the HT interface.
This example show how to build a flow to compute the band structure and the DOS of silicon.
Input files are automatically generated with factory functions designed for automatic calculations.
"""
from __future__ import print_function, division, unicode_literals, absolute_import
@ -17,7 +18,7 @@ from abipy import abilab
def build_flow(options):
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_","flow_")
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
# Initialize structure and pseudos.
structure = abilab.Structure.from_file(abidata.cif_file("si.cif"))
@ -26,11 +27,14 @@ def build_flow(options):
# Initialize the flow.
flow = flowtk.Flow(workdir=options.workdir, manager=options.manager)
# Use ebands_input factory function to build inputs.
multi = abilab.ebands_input(structure, pseudos, kppa=40, nscf_nband=6, ndivsm=10, ecut=6)
work = flowtk.BandStructureWork(scf_input=multi[0], nscf_input=multi[1])
# Use the ebands_input factory function to build a MultiDataset.
# keyword args are optional (default values are given or computed automatically, see docs).
multi = abilab.ebands_input(structure, pseudos, kppa=40, dos_kppa=80,
nscf_nband=6, ndivsm=10, ecut=6, spin_mode="unpolarized")
work = flowtk.BandStructureWork(scf_input=multi[0], nscf_input=multi[1], dos_inputs=multi[2])
flow.register_work(work)
return flow
@ -40,7 +44,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -55,3 +59,42 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_hs_si_ebands.py -s
#
# Use:
#
# abirun.py flow_h_si_ebands deps
#
# To print the connections among the tasks.
#
# .. code-block:: bash
#
# <ScfTask, node_id=241095, workdir=flow_ht_si_ebands/w0/t0>
#
# <NscfTask, node_id=241096, workdir=flow_ht_si_ebands/w0/t1>
# +--<ScfTask, node_id=241095, workdir=flow_ht_si_ebands/w0/t0>
#
# <NscfTask, node_id=241097, workdir=flow_ht_si_ebands/w0/t2>
# +--<ScfTask, node_id=241095, workdir=flow_ht_si_ebands/w0/t0>
#
# ``w0/t1`` is a band structure calculation along a k-path
# while ``w0/t2`` produced KS eigenvalues in the IBZ for the DOS.
#
# You can use ipython to plot the electronic bands with the DOS:
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: from abipy import abilab
# In [3]: gsr_kpath = abilab.abiopen("flow_ht_si_ebands/w0/t1/outdata/out_GSR.nc")
# In [4]: gsr_kmesh = abilab.abiopen("flow_ht_si_ebands/w0/t2/outdata/out_GSR.nc")
# In [5]: gsr_kpath.ebands.plot_with_edos(gsr_kmesh.ebands.get_edos())
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_ht_si_ebands.png?raw=true
# :alt: Band structure of Si in the IBZ and along a k-path
#

160
abipy/examples/flows/run_ht_si_g0w0ppm.py
View File

@ -21,34 +21,30 @@ def build_flow(options):
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_","flow_")
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
# Initialize the flow.
flow = flowtk.Flow(options.workdir, manager=options.manager)
scf_kppa = 10
nscf_nband = 10
#nscf_ngkpt = [4,4,4]
#nscf_shiftk = [0.0, 0.0, 0.0]
ecut, ecuteps, ecutsigx = 4, 2, 3
scf_kppa = 120
nscf_nband = 40
ecut, ecuteps, ecutsigx = 6, 2, 4
#scr_nband = 50
#sigma_nband = 50
multi = abilab.g0w0_with_ppmodel_inputs(
structure, pseudos,
scf_kppa, nscf_nband, ecuteps, ecutsigx,
ecut=ecut, pawecutdg=None,
structure, pseudos, scf_kppa, nscf_nband, ecuteps, ecutsigx,
ecut=ecut, shifts=(0, 0, 0), # By default the k-mesh is shifted! TODO: Change default?
accuracy="normal", spin_mode="unpolarized", smearing=None,
#ppmodel="godby", charge=0.0, scf_algorithm=None, inclvkb=2, scr_nband=None,
#sigma_nband=None, gw_qprange=1):
)
#multi.set_vars(paral_kgb=1)
scf_input, nscf_input, scr_input, sigma_input = multi.split_datasets()
work = flowtk.G0W0Work(scf_input, nscf_input, scr_input, sigma_input)
flow.register_work(work)
return flow
@ -58,7 +54,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -73,3 +69,143 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_ht_si_g0w0ppm.py -s
#
# The last task ``w0_t3`` is the SigmaTask who has produced a netcdf file with the GW results
# as you can we with the command:
#
# abirun.py flow_ht_si_g0w0ppm/ listext SIGRES
#
# Let's print the QP results with:
#
# abiopen.py flow_ht_si_g0w0ppm/w0/t3/outdata/out_SIGRES.nc -p
#
# .. code-block:: bash
#
# ================================= File Info =================================
# Name: out_SIGRES.nc
# Directory: /Users/gmatteo/git_repos/abipy/abipy/examples/flows/flow_ht_si_g0w0ppm/w0/t3/outdata
# Size: 254.07 kb
# Access Time: Sat Dec 9 19:16:05 2017
# Modification Time: Sat Dec 9 15:16:15 2017
# Change Time: Sat Dec 9 15:16:15 2017
#
# ================================= Structure =================================
# Full Formula (Si2)
# Reduced Formula: Si
# abc : 3.866975 3.866975 3.866975
# angles: 60.000000 60.000000 60.000000
# Sites (2)
# # SP a b c
# --- ---- ---- ---- ----
# 0 Si 0 0 0
# 1 Si 0.25 0.25 0.25
#
# Abinit Spacegroup: spgid: 0, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True
#
# ============================== Kohn-Sham bands ==============================
# Number of electrons: 8.0, Fermi level: 5.963 [eV]
# nsppol: 1, nkpt: 8, mband: 40, nspinor: 1, nspden: 1
# smearing scheme: none, tsmear_eV: 0.272, occopt: 1
# Direct gap:
# Energy: 2.512 [eV]
# Initial state: spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.016, band=3, eig=5.640, occ=2.000
# Final state: spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.016, band=4, eig=8.152, occ=0.000
# Fundamental gap:
# Energy: 0.646 [eV]
# Initial state: spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.016, band=3, eig=5.640, occ=2.000
# Final state: spin=0, kpt=[+0.500, +0.500, +0.000], weight: 0.047, band=4, eig=6.286, occ=0.000
# Bandwidth: 11.867 [eV]
# Valence minimum located at:
# spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.016, band=0, eig=-6.227, occ=2.000
# Valence maximum located at:
# spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.016, band=3, eig=5.640, occ=2.000
#
# ============================ QP direct gaps in eV ============================
# QP_dirgap: 3.000 for K-point: [+0.000, +0.000, +0.000], spin: 0
# QP_dirgap: 3.087 for K-point: [+0.250, +0.000, +0.000], spin: 0
# QP_dirgap: 3.068 for K-point: [+0.500, +0.000, +0.000], spin: 0
# QP_dirgap: 3.421 for K-point: [+0.250, +0.250, +0.000], spin: 0
# QP_dirgap: 4.165 for K-point: [+0.500, +0.250, +0.000], spin: 0
# QP_dirgap: 4.206 for K-point: [-0.250, +0.250, +0.000], spin: 0
# QP_dirgap: 4.003 for K-point: [+0.500, +0.500, +0.000], spin: 0
# QP_dirgap: 8.683 for K-point: [-0.250, +0.500, +0.250], spin: 0
#
# ============== QP results for each k-point and spin (All in eV) ==============
# K-point: [+0.000, +0.000, +0.000], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 1 1 5.640 5.547 5.520 -11.156 -12.870 1.594 0.0 0.776
# 2 2 5.640 5.547 5.521 -11.156 -12.870 1.594 0.0 0.776
# 3 3 5.640 5.547 5.521 -11.156 -12.870 1.594 0.0 0.776
# 4 4 8.152 8.547 8.662 -9.969 -5.622 -3.837 0.0 0.774
# 5 5 8.152 8.547 8.662 -9.969 -5.622 -3.837 0.0 0.774
# 6 6 8.152 8.547 8.662 -9.969 -5.622 -3.837 0.0 0.774
#
# K-point: [+0.250, +0.000, +0.000], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 2 2 4.870 4.751 4.716 -10.936 -12.831 1.741 0.0 0.772
# 3 3 4.870 4.751 4.716 -10.936 -12.831 1.741 0.0 0.772
# 4 4 7.449 7.838 7.949 -10.031 -5.725 -3.806 0.0 0.778
# 5 5 9.108 9.534 9.660 -10.013 -5.381 -4.079 0.0 0.771
# 6 6 9.108 9.534 9.660 -10.013 -5.381 -4.079 0.0 0.771
#
# K-point: [+0.500, +0.000, +0.000], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 2 2 4.431 4.282 4.238 -10.908 -13.077 1.976 0.0 0.769
# 3 3 4.431 4.282 4.238 -10.908 -13.077 1.976 0.0 0.769
# 4 4 6.972 7.350 7.456 -10.016 -5.844 -3.687 0.0 0.781
# 5 5 8.982 9.410 9.534 -9.621 -4.942 -4.127 0.0 0.775
# 6 6 8.982 9.410 9.534 -9.621 -4.942 -4.127 0.0 0.775
#
# K-point: [+0.250, +0.250, +0.000], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 2 2 3.743 3.571 3.518 -10.603 -12.883 2.056 0.0 0.764
# 3 3 3.743 3.571 3.518 -10.603 -12.883 2.056 0.0 0.764
# 4 4 6.692 6.993 7.076 -9.352 -5.373 -3.595 0.0 0.783
# 5 5 8.669 9.035 9.139 -9.203 -4.548 -4.184 0.0 0.779
#
# K-point: [+0.500, +0.250, +0.000], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 2 2 2.099 1.893 1.825 -10.130 -12.839 2.435 0.0 0.752
# 3 3 3.427 3.242 3.184 -10.596 -13.077 2.238 0.0 0.761
# 4 4 7.086 7.408 7.496 -9.224 -5.065 -3.748 0.0 0.784
# 5 5 10.015 10.397 10.510 -9.594 -4.730 -4.369 0.0 0.771
#
# K-point: [-0.250, +0.250, +0.000], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 2 2 1.887 1.667 1.594 -9.982 -12.769 2.493 0.0 0.751
# 3 3 4.302 4.157 4.113 -10.820 -12.927 1.917 0.0 0.769
# 4 4 7.983 8.363 8.470 -9.567 -5.113 -3.966 0.0 0.779
# 5 5 10.350 10.800 10.942 -10.046 -4.859 -4.595 0.0 0.760
#
# K-point: [+0.500, +0.500, +0.000], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 2 2 2.784 2.567 2.497 -10.480 -13.270 2.503 0.0 0.755
# 3 3 2.784 2.567 2.497 -10.480 -13.270 2.503 0.0 0.755
# 4 4 6.286 6.571 6.647 -9.003 -5.038 -3.603 0.0 0.787
# 5 5 6.286 6.571 6.647 -9.003 -5.038 -3.603 0.0 0.787
#
# K-point: [-0.250, +0.500, +0.250], spin: 0
# band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
# 2 2 1.787 1.557 1.481 -9.917 -12.751 2.527 0.0 0.749
# 3 3 1.787 1.557 1.481 -9.917 -12.751 2.527 0.0 0.749
# 4 4 9.875 10.241 10.347 -9.523 -4.876 -4.176 0.0 0.775
# 5 5 9.875 10.241 10.347 -9.523 -4.876 -4.176 0.0 0.775
#
# we can also plot the QP data as function of the initial KS energy with ipython:
#
# abiopen.py flow_ht_si_g0w0ppm/w0/t3/outdata/out_SIGRES.nc
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: abifile.plot_qps_vs_e0()
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_ht_si_g0w0_ppm.png?raw=true
# :alt: QP results in Si plotted vs the KS energy e0.
#

71
abipy/examples/flows/run_ldaus.py
View File

@ -3,7 +3,8 @@ r"""
Flow for LDA+U calculations
===========================
LDA+U band structure of NiO for several values of U-J.
This example shows how to compute the LDA+U band structure of NiO
with PAW for several values of U-J.
"""
from __future__ import print_function, division, unicode_literals, absolute_import
@ -58,7 +59,7 @@ def make_scf_nscf_dos_inputs(structure, pseudos, luj_params, paral_kgb=1):
multi[1].set_kpath(ndivsm=6)
multi[1].set_vars(tolwfr=1e-10)
# Dos calculation.
# DOS calculation.
multi[2].set_vars(
iscf=-3, # NSCF calculation
ngkpt=structure.calc_ngkpt(nksmall=8),
@ -112,7 +113,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -127,3 +128,67 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
# Run the script with:
#
# run_ldaus.py -s
#
# then use:
#
# abirun.py flow_ldaus/ ebands
#
# to analyze all GSR files produced by flow.
#
# .. code-block:: bash
#
# KS electronic bands:
# nsppol nspinor nspden nkpt nband nelect fermie formula natom \
# w0_t0 1 1 2 3 40 48.0 6.084 Ni2 O2 4
# w0_t1 1 1 2 129 40 48.0 6.084 Ni2 O2 4
# w0_t2 1 1 2 213 40 48.0 6.169 Ni2 O2 4
# w1_t0 1 1 2 3 40 48.0 6.855 Ni2 O2 4
# w1_t1 1 1 2 129 40 48.0 6.855 Ni2 O2 4
# w1_t2 1 1 2 213 40 48.0 6.432 Ni2 O2 4
#
# angle0 angle1 angle2 a b c volume abispg_num \
# w0_t0 60.0 60.0 60.0 2.964 2.964 5.927 36.809 166
# w0_t1 60.0 60.0 60.0 2.964 2.964 5.927 36.809 166
# w0_t2 60.0 60.0 60.0 2.964 2.964 5.927 36.809 166
# w1_t0 60.0 60.0 60.0 2.964 2.964 5.927 36.809 166
# w1_t1 60.0 60.0 60.0 2.964 2.964 5.927 36.809 166
# w1_t2 60.0 60.0 60.0 2.964 2.964 5.927 36.809 166
#
# scheme occopt tsmear_ev bandwidth_spin0 fundgap_spin0 \
# w0_t0 gaussian 7 0.408 99.343 3.324
# w0_t1 gaussian 7 0.408 99.680 3.049
# w0_t2 gaussian 7 0.408 99.838 2.560
# w1_t0 gaussian 7 0.408 99.318 4.626
# w1_t1 gaussian 7 0.408 99.628 3.352
# w1_t2 gaussian 7 0.408 99.826 3.155
#
# dirgap_spin0 task_class ncfile node_id \
# w0_t0 3.351 ScfTask flow_ldaus/w0/t0/outdata/out_GSR.nc 241313
# w0_t1 3.352 NscfTask flow_ldaus/w0/t1/outdata/out_GSR.nc 241314
# w0_t2 3.041 NscfTask flow_ldaus/w0/t2/outdata/out_GSR.nc 241315
# w1_t0 4.626 ScfTask flow_ldaus/w1/t0/outdata/out_GSR.nc 241317
# w1_t1 3.928 NscfTask flow_ldaus/w1/t1/outdata/out_GSR.nc 241318
# w1_t2 3.983 NscfTask flow_ldaus/w1/t2/outdata/out_GSR.nc 241319
#
# status
# w0_t0 Completed
# w0_t1 Completed
# w0_t2 Completed
# w1_t0 Completed
# w1_t1 Completed
# w1_t2 Completed
#
# The second task in each work (w*_t1) is a NSCF run along the high symmetry path.
# and we can compare the results of these two task by specifying the node identifiers:
#
# abirun.py flow_ldaus/ ebands -p --nids=241314,241318
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_ldaus.png?raw=true
# :alt: Band structure of Si in the IBZ and along a k-path
#

97
abipy/examples/flows/run_mgb2_edoses.py
View File

@ -1,9 +1,10 @@
#!/usr/bin/env python
r"""
Flow for Bands + DOS
====================
Flow for convergence studies of e-DOS wrt ngkpt
===============================================
Band structure and the electron DOS of MgB2 with different k-point samplings.
This examples shows how to build a Flow to compute the
band structure and the electron DOS of MgB2 with different k-point samplings.
"""
from __future__ import print_function, division, unicode_literals, absolute_import
@ -20,15 +21,14 @@ def make_scf_nscf_inputs(structure, pseudos, paral_kgb=1):
multi = abilab.MultiDataset(structure, pseudos=pseudos, ndtset=5)
# Global variables
global_vars = dict(ecut=10,
nband=11,
timopt=-1,
occopt=4, # Marzari smearing
tsmear=0.03,
paral_kgb=paral_kgb,
)
multi.set_vars(global_vars)
multi.set_vars(
ecut=10,
nband=11,
timopt=-1,
occopt=4, # Marzari smearing
tsmear=0.03,
paral_kgb=paral_kgb,
)
# Dataset 1 (GS run)
multi[0].set_kmesh(ngkpt=[8,8,8], shiftk=structure.calc_shiftk())
@ -68,7 +68,6 @@ def build_flow(options):
inputs = make_scf_nscf_inputs(structure, pseudos)
scf_input, nscf_input, dos_inputs = inputs[0], inputs[1], inputs[2:]
#print(scf_input.pseudos)
return flowtk.bandstructure_flow(options.workdir, scf_input, nscf_input,
dos_inputs=dos_inputs, manager=options.manager)
@ -80,7 +79,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -95,3 +94,73 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_mgb2_edoses.py -s
#
# then use:
#
# abirun.py flow_mgb2_edoes ebands
#
# to get info about the electronic properties:
#
# .. code-block:: shell
#
# KS electronic bands:
# nsppol nspinor nspden nkpt nband nelect fermie formula natom \
# w0_t0 1 1 1 40 11 8.0 7.615 Mg1 B2 3
# w0_t1 1 1 1 97 11 8.0 7.615 Mg1 B2 3
# w0_t2 1 1 1 15 11 8.0 7.701 Mg1 B2 3
# w0_t3 1 1 1 80 11 8.0 7.629 Mg1 B2 3
# w0_t4 1 1 1 432 11 8.0 7.626 Mg1 B2 3
#
# angle0 angle1 angle2 a b c volume abispg_num \
# w0_t0 90.0 90.0 120.0 3.086 3.086 3.523 29.056 191
# w0_t1 90.0 90.0 120.0 3.086 3.086 3.523 29.056 191
# w0_t2 90.0 90.0 120.0 3.086 3.086 3.523 29.056 191
# w0_t3 90.0 90.0 120.0 3.086 3.086 3.523 29.056 191
# w0_t4 90.0 90.0 120.0 3.086 3.086 3.523 29.056 191
#
# scheme occopt tsmear_ev \
# w0_t0 cold smearing of N. Marzari with minimization ... 4 0.816
# w0_t1 cold smearing of N. Marzari with minimization ... 4 0.816
# w0_t2 cold smearing of N. Marzari with minimization ... 4 0.816
# w0_t3 cold smearing of N. Marzari with minimization ... 4 0.816
# w0_t4 cold smearing of N. Marzari with minimization ... 4 0.816
#
# bandwidth_spin0 fundgap_spin0 dirgap_spin0 task_class \
# w0_t0 12.452 0.031 0.609 ScfTask
# w0_t1 12.441 0.077 0.399 NscfTask
# w0_t2 12.368 0.415 1.680 NscfTask
# w0_t3 12.510 0.069 0.390 NscfTask
# w0_t4 12.506 0.033 0.283 NscfTask
#
# ncfile node_id status
# w0_t0 flow_mgb2_edoses/w0/t0/outdata/out_GSR.nc 241032 Completed
# w0_t1 flow_mgb2_edoses/w0/t1/outdata/out_GSR.nc 241033 Completed
# w0_t2 flow_mgb2_edoses/w0/t2/outdata/out_GSR.nc 241034 Completed
# w0_t3 flow_mgb2_edoses/w0/t3/outdata/out_GSR.nc 241035 Completed
# w0_t4 flow_mgb2_edoses/w0/t4/outdata/out_GSR.nc 241036 Completed
#
# Our main goal is to analyze the convergence of the DOS wrt to the k-point sampling.
# As we know that ``w0_t2``, `w0_t3`` and ``w0_t4`` are DOS calculations, we can
# build a GSR robot for these tasks with:
#
# abirun.py flow_mgb2_edoses/ ebands -nids=241034,241035,241036
#
# then, inside the ipython shell, one can use:
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: robot.combiplot_edos()
#
# to plot the electronic DOS obtained with different number of k-points in the IBZ.
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_mgb2_edoses.png?raw=true
# :alt: Convergence of electronic DOS in MgB2 wrt k-points.
#

2
abipy/examples/flows/run_nonlinear.py
View File

@ -68,7 +68,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=False, tight_layout=True)

143
abipy/examples/flows/run_optic.py
View File

@ -3,7 +3,9 @@ r"""
Optic Flow
==========
Optical spectra with Optic.
This example shows how to create a flow to compute optical spectra with Optic
(independent particle approximation, no local field effects) and perform
a convergence study with respect to the k-point sampling.
"""
from __future__ import print_function, division, unicode_literals, absolute_import
@ -15,116 +17,74 @@ import abipy.flowtk as flowtk
def build_flow(options, paral_kgb=0):
"""
Build flow for the calculation of optical properties with optic + band structure
along high-symmetry k-path. DDK are computed with 3 k-meshes of increasing density
to monitor the convergece of the spectra.
"""
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
multi = abilab.MultiDataset(structure=abidata.structure_from_ucell("GaAs"),
pseudos=abidata.pseudos("31ga.pspnc", "33as.pspnc"), ndtset=5)
pseudos=abidata.pseudos("31ga.pspnc", "33as.pspnc"), ndtset=2)
# Global variables
kmesh = dict(ngkpt=[4, 4, 4],
nshiftk=4,
shiftk=[[0.5, 0.5, 0.5],
[0.5, 0.0, 0.0],
[0.0, 0.5, 0.0],
[0.0, 0.0, 0.5]]
)
# Usa same shifts in all tasks.
shiftk= [[0.5, 0.5, 0.5],
[0.5, 0.0, 0.0],
[0.0, 0.5, 0.0],
[0.0, 0.0, 0.5]]
global_vars = dict(ecut=2, paral_kgb=paral_kgb)
global_vars.update(kmesh)
multi.set_vars(global_vars)
# Global variables.
multi.set_vars(ecut=2, paral_kgb=paral_kgb)
# Dataset 1 (GS run)
multi[0].set_vars(
tolvrs=1e-6,
nband=4,
)
multi[0].set_vars(tolvrs=1e-8, nband=4)
multi[0].set_kmesh(ngkpt=[4, 4, 4], shiftk=shiftk)
# NSCF run with large number of bands, and points in the the full BZ
multi[1].set_vars(
iscf=-2,
nband=20,
nstep=25,
kptopt=1,
tolwfr=1.e-9,
#kptopt=3,
)
# Fourth dataset : ddk response function along axis 1
# Fifth dataset : ddk response function along axis 2
# Sixth dataset : ddk response function along axis 3
for idir in range(3):
rfdir = 3 * [0]
rfdir[idir] = 1
multi[2+idir].set_vars(
iscf=-3,
nband=20,
nstep=1,
nline=0,
prtwf=3,
kptopt=3,
nqpt=1,
qpt=[0.0, 0.0, 0.0],
rfdir=rfdir,
rfelfd=2,
tolwfr=1.e-9,
)
scf_inp, nscf_inp, ddk1, ddk2, ddk3 = multi.split_datasets()
# NSCF run on k-path with large number of bands
multi[1].set_vars(iscf=-2, nband=20, tolwfr=1.e-9)
multi[1].set_kpath(ndivsm=10)
# Initialize the flow.
flow = flowtk.Flow(options.workdir, manager=options.manager)
# GS to get the density + NSCF along the path.
scf_inp, nscf_inp = multi.split_datasets()
bands_work = flowtk.BandStructureWork(scf_inp, nscf_inp)
flow.register_work(bands_work)
ddk_work = flowtk.Work()
for inp in [ddk1, ddk2, ddk3]:
ddk_work.register_ddk_task(inp, deps={bands_work.nscf_task: "WFK"})
flow.register_work(ddk_work)
# Optic does not support MPI with ncpus > 1.
# Build OpticInput used to compute optical properties.
optic_input = abilab.OpticInput(
broadening=0.002, # Value of the smearing factor, in Hartree
domega=0.0003, # Frequency mesh.
maxomega=0.3,
scissor=0.000, # Scissor shift if needed, in Hartree
tolerance=0.002, # Tolerance on closeness of singularities (in Hartree)
num_lin_comp=1, # Number of components of linear optic tensor to be computed
lin_comp=11, # Linear coefficients to be computed (x=1, y=2, z=3)
num_lin_comp=2, # Number of components of linear optic tensor to be computed
lin_comp=(11, 33), # Linear coefficients to be computed (x=1, y=2, z=3)
num_nonlin_comp=2, # Number of components of nonlinear optic tensor to be computed
nonlin_comp=(123, 222), # Non-linear coefficients to be computed
)
# TODO
# Check is the order of the 1WF files is relevant. Can we use DDK files ordered
# in an arbitrary way or do we have to pass (x,y,z)?
optic_task = flowtk.OpticTask(optic_input, nscf_node=bands_work.nscf_task, ddk_nodes=ddk_work)
flow.register_task(optic_task)
# ddk_nband is fixed here, in principle it depends on nelect and the frequency range in chi(w).
ddk_nband = 20
return flow
# Perform converge study wrt ngkpt (shiftk is constant).
ngkpt_convergence = [[4, 4, 4], [8, 8, 8], [16, 16, 16]]
from abipy.flowtk.dfpt_works import NscfDdksWork
for ddk_ngkpt in ngkpt_convergence:
# Build work for NSCF from DEN produced by the first GS task + 3 DDKs.
# All tasks use more bands and a denser k-mesh defined by ddk_ngkpt.
ddks_work = NscfDdksWork.from_scf_task(bands_work[0], ddk_ngkpt, shiftk, ddk_nband)
flow.register_work(ddks_work)
def optic_flow_from_files():
# Optic does not support MPI with ncpus > 1.
manager = flowtk.TaskManager.from_user_config()
manager.set_mpi_procs(1)
# Build optic task to compute chi with this value of ddk_ngkpt.
optic_task = flowtk.OpticTask(optic_input, nscf_node=ddks_work.task_with_ks_energies,
ddk_nodes=ddks_work.ddk_tasks, use_ddknc=False)
ddks_work.register_task(optic_task)
flow = flowtk.Flow(workdir="OPTIC_FROM_FILE", manager=manager)
ddk_nodes = [
"/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_1/task_0/outdata/out_1WF",
"/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_1/task_1/outdata/out_1WF",
"/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_1/task_2/outdata/out_1WF",
]
nscf_node = "/Users/gmatteo/Coding/abipy/abipy/data/runs/OPTIC/work_0/task_1/outdata/out_WFK"
optic_task = flowtk.OpticTask(optic_input, nscf_node=nscf_node, ddk_nodes=ddk_nodes)
flow.register_task(optic_task)
return flow
@ -134,7 +94,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=False, tight_layout=True)
@flowtk.flow_main
@ -149,3 +109,24 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_optic.py -s
#
# then use:
#
# abirun.py flow_optic robot optic
#
# to create a robot for OPTIC.nc files. Then inside the ipytho shell type:
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: robot.plot_linopt_convergence()
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_optic.png?raw=true
# :alt: Convergence of (linear) optical spectra wrt k-points.
#

108
abipy/examples/flows/run_optic_new.py
View File

@ -1,108 +0,0 @@
#!/usr/bin/env python
r"""
Optic Flow
==========
Optical spectra with Optic.
"""
from __future__ import print_function, division, unicode_literals, absolute_import
import sys
import os
import abipy.data as abidata
import abipy.abilab as abilab
import abipy.flowtk as flowtk
def build_flow(options, paral_kgb=0):
"""
Build flow for the calculation of optical properties with optic + band structure
along high-symmetry k-path. DDK are computed with 3 k-meshes of increasing density
to monitor the convergece of the spectra.
"""
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
multi = abilab.MultiDataset(structure=abidata.structure_from_ucell("GaAs"),
pseudos=abidata.pseudos("31ga.pspnc", "33as.pspnc"), ndtset=2)
# Usa same shifts in all tasks.
shiftk= [[0.5, 0.5, 0.5],
[0.5, 0.0, 0.0],
[0.0, 0.5, 0.0],
[0.0, 0.0, 0.5]]
# Global variables.
multi.set_vars(ecut=2, paral_kgb=paral_kgb)
# Dataset 1 (GS run)
multi[0].set_vars(tolvrs=1e-8, nband=4)
multi[0].set_kmesh(ngkpt=[4, 4, 4], shiftk=shiftk)
# NSCF run on k-path with large number of bands
multi[1].set_vars(iscf=-2, nband=20, tolwfr=1.e-9)
multi[1].set_kpath(ndivsm=10)
# Initialize the flow.
flow = flowtk.Flow(options.workdir, manager=options.manager)
# GS to get the density + NSCF along the path.
scf_inp, nscf_inp = multi.split_datasets()
bands_work = flowtk.BandStructureWork(scf_inp, nscf_inp)
flow.register_work(bands_work)
# Build OpticInput used to compute optical properties.
optic_input = abilab.OpticInput(
broadening=0.002, # Value of the smearing factor, in Hartree
domega=0.0003, # Frequency mesh.
maxomega=0.3,
scissor=0.000, # Scissor shift if needed, in Hartree
tolerance=0.002, # Tolerance on closeness of singularities (in Hartree)
num_lin_comp=2, # Number of components of linear optic tensor to be computed
lin_comp=(11, 33), # Linear coefficients to be computed (x=1, y=2, z=3)
num_nonlin_comp=2, # Number of components of nonlinear optic tensor to be computed
nonlin_comp=(123, 222), # Non-linear coefficients to be computed
)
# ddk_nband is fixed here, in principle it depends on nelect and the frequency range in chi(w).
ddk_nband = 20
# Perform converge study wrt ngkpt (shiftk is constant).
ngkpt_convergence = [[4, 4, 4], [8, 8, 8], [16, 16, 16]]
for ddk_ngkpt in ngkpt_convergence:
# Build work for NSCF from DEN produced by the first GS task + 3 DDKs.
# All tasks use more bands and a denser k-mesh defined by ddk_ngkpt.
ddks_work = flowtk.NscfDdksWork.from_scf_task(bands_work[0], ddk_ngkpt, shiftk, ddk_nband)
flow.register_work(ddks_work)
# Build optic task to compute chi with this value of ddk_ngkpt.
optic_task = flowtk.OpticTask(optic_input, nscf_node=ddks_work.task_with_ks_energies,
ddk_nodes=ddks_work.ddk_tasks, use_ddknc=False)
ddks_work.register_task(optic_task)
return flow
# This block generates the thumbnails in the Abipy gallery.
# You can safely REMOVE this part if you are using this script for production runs.
if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
@flowtk.flow_main
def main(options):
"""
This is our main function that will be invoked by the script.
flow_main is a decorator implementing the command line interface.
Command line args are stored in `options`.
"""
return build_flow(options)
if __name__ == "__main__":
sys.exit(main())

73
abipy/examples/flows/run_phfrozen_ebands.py
View File

@ -85,7 +85,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -100,3 +100,74 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_phfrozen_ebands -s
#
# then use:
#
# abirun.py flow_phfrozen_ebands/ structures -v
#
# to analyze the input/output structures including the atomic positions:
#
# .. code-block:: bash
#
# Lattice parameters:
# formula natom angle0 angle1 angle2 a b c volume \
# w0_t0_in Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w0_t0_out Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w0_t1_in Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w1_t0_in Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w1_t0_out Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w1_t1_in Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w2_t0_in Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w2_t0_out Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
# w2_t1_in Si2 2 60.0 60.0 60.0 3.867 3.867 3.867 40.888
#
# abispg_num P [GPa] Max|F| eV/ang task_class status
# w0_t0_in None NaN NaN ScfTask Completed
# w0_t0_out 12 -3.638 3.180e+00 ScfTask Completed
# w0_t1_in None NaN NaN NscfTask Completed
# w1_t0_in None NaN NaN ScfTask Completed
# w1_t0_out 227 -5.212 7.430e-27 ScfTask Completed
# w1_t1_in None NaN NaN NscfTask Completed
# w2_t0_in None NaN NaN ScfTask Completed
# w2_t0_out 12 -4.192 2.095e+00 ScfTask Completed
# w2_t1_in None NaN NaN NscfTask Completed
#
# Atomic positions (columns give the site index):
# 0 \
# w0_t0_in (Si, +0.970139 +0.000000 +0.014930)
# w0_t0_out (Si, +0.970139 +0.000000 +0.014930)
# w0_t1_in (Si, +0.970139 +0.000000 +0.014930)
# w1_t0_in (Si, +0.000000 +0.000000 +0.000000)
# w1_t0_out (Si, +0.000000 +0.000000 +0.000000)
# w1_t1_in (Si, +0.000000 +0.000000 +0.000000)
# w2_t0_in (Si, +0.029861 +0.000000 +0.985070)
# w2_t0_out (Si, +0.029861 +0.000000 +0.985070)
# w2_t1_in (Si, +0.029861 +0.000000 +0.985070)
#
# 1 status
# w0_t0_in (Si, +0.279861 +0.250000 +0.235070) Completed
# w0_t0_out (Si, +0.279861 +0.250000 +0.235070) Completed
# w0_t1_in (Si, +0.279861 +0.250000 +0.235070) Completed
# w1_t0_in (Si, +0.250000 +0.250000 +0.250000) Completed
# w1_t0_out (Si, +0.250000 +0.250000 +0.250000) Completed
# w1_t1_in (Si, +0.250000 +0.250000 +0.250000) Completed
# w2_t0_in (Si, +0.220139 +0.250000 +0.264930) Completed
# w2_t0_out (Si, +0.220139 +0.250000 +0.264930) Completed
# w2_t1_in (Si, +0.220139 +0.250000 +0.264930) Completed
#
# Finally, we can plot the electronic bands with the command:
#
# abirun.py flow_phfrozen_ebands ebands -p -t NscfTask
#
# to select only the band structures produced by the NscfTask.
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_phfrozen_ebands.png?raw=true
# :alt: Band structures of Si computed for different displacement amplitudes.
#

47
abipy/examples/flows/run_phonons.py
View File

@ -3,7 +3,11 @@ r"""
Flow for phonons with DFPT
==========================
Phonon band structure of AlAs.
This example shows how to compute the phonon band structure of AlAs with AbiPy flows.
Symmetries are taken into account: only q-points in the IBZ are generated and
for each q-point only the independent atomic perturbations are computed.
The final results (out_DDB, out_DVDB) will be produced automatically at the end of the run
and saved in ``flow_phonons/w1/outdata/``.
"""
from __future__ import division, print_function, unicode_literals, absolute_import
@ -37,7 +41,6 @@ def make_scf_input(paral_kgb=0):
paral_kgb=paral_kgb,
tolvrs=1.0e-10,
ixc=1,
nstep=25,
diemac=9.0,
iomode=3,
)
@ -79,7 +82,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=False, tight_layout=True)
@flowtk.flow_main
@ -94,3 +97,41 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_phonons.py -s
#
# then use:
#
# abirun.py flow_phonons history
#
# to get the list of actions perfomed by AbiPy to complete the flow.
# Note how the ``PhononWork`` has merged all the partial DDB files produced by the PhononTasks
#
# .. code-block:: bash
#
# ===================================================================================================================================
# ====================================== <PhononWork, node_id=241274, workdir=flow_phonons/w1> ======================================
# ===================================================================================================================================
# [Thu Dec 7 22:55:02 2017] Finalized set to True
# [Thu Dec 7 22:55:02 2017] Will call mrgddb to merge [ .... ]
#
# Now open the final DDB file with:
#
# abiopen.py flow_phonons/w1/outdata/out_DDB
#
# and invoke anaddb to compute the phonon band structure and the phonon DOS with:
#
# .. code-block:: ipython
#
# In [1]: phbst_file, phdos_file = abifile.anaget_phbst_and_phdos_files()
# In [2]: %matplotlib
# In [3]: phbst_file.plot_phbands()
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_phonons.png?raw=true
# :alt: Phonon band structure of AlAs.
#

44
abipy/examples/flows/run_phonopy_si.py
View File

@ -3,7 +3,12 @@ r"""
Phonopy + Abinit Flow
=====================
Compute phonon frequencies with phonopy (supercells and finite-difference method).
This example shows how to compute phonon frequencies with phonopy (supercells and finite-difference method).
This approach could be useful to obtain vibrational properties with XC functionals for which DFPT is not yet implemented.
.. warning:
This example requires the `phonopy package <http://atztogo.github.io/phonopy/examples.html>`_
"""
from __future__ import print_function, division, unicode_literals, absolute_import
@ -31,6 +36,7 @@ def build_flow(options):
# Build input for GS calculation.
gsinp = abilab.AbinitInput(structure, pseudos)
gsinp.set_vars(ecut=4, nband=4, toldff=1.e-6)
# This gives ngkpt = 4x4x4 with 4 shifts for the initial unit cell.
# The k-point sampling will be rescaled when we build the supercell in PhonopyWork.
gsinp.set_autokmesh(nksmall=4)
@ -51,7 +57,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -66,3 +72,37 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_phonopy_si.py -s
#
# the output results are produced in ``flow_phonopy_si/w0/outdata/``
# Follow the instructions in the README file:
#
# .. code-block:: md
#
# To plot bands, use:
# phonopy -p band.conf
#
# To plot phonon dos, use:
# phonopy -p dos.conf
#
# To plot bands and dos, use:
# phonopy -p band-dos.conf
#
# See also:
# http://atztogo.github.io/phonopy/examples.html
# http://atztogo.github.io/phonopy/setting-tags.html#setting-tags
#
# The command:
#
# phonopy -p band-dos.conf
#
# will produce:
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_phonopy_si.png?raw=true
# :alt: Phonon Band structure computed with phonopy.
#

44
abipy/examples/flows/run_raman_bse.py
View File

@ -69,25 +69,26 @@ def raman_work(structure, pseudos, shiftk, paral_kgb=1):
# GS run
scf_inp = abilab.AbinitInput(structure, pseudos=pseudos)
scf_inp.set_vars(global_vars)
scf_inp.set_kmesh(ngkpt=[2,2,2], shiftk=shiftk)
scf_inp.set_kmesh(ngkpt=[2, 2, 2], shiftk=shiftk)
scf_inp["tolvrs"] = 1e-6
# NSCF run
nscf_inp = abilab.AbinitInput(structure, pseudos=pseudos)
nscf_inp.set_vars(global_vars)
nscf_inp.set_kmesh(ngkpt=[2,2,2], shiftk=shiftk)
nscf_inp.set_kmesh(ngkpt=[2, 2, 2], shiftk=shiftk)
nscf_inp.set_vars(tolwfr=1e-8,
nband=12,
nbdbuf=4,
iscf=-2,
)
nscf_inp.set_vars(
tolwfr=1e-8,
nband=12,
nbdbuf=4,
iscf=-2,
)
# BSE run with Model dielectric function and Haydock (only resonant + W + v)
# Note that SCR file is not needed here
bse_inp = abilab.AbinitInput(structure, pseudos=pseudos)
bse_inp.set_vars(global_vars)
bse_inp.set_kmesh(ngkpt=[2,2,2], shiftk=shiftk)
bse_inp.set_kmesh(ngkpt=[2, 2, 2], shiftk=shiftk)
bse_inp.set_vars(
optdriver=99,
@ -118,7 +119,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main
@ -133,3 +134,28 @@ def main(options):
if __name__ == "__main__":
sys.exit(main())
############################################################################
#
# Run the script with:
#
# run_raman_bse.py -s
#
# then use:
#
# abirun.py flow_raman_bse robot MDF
#
# to analyze all the macroscopic dielectric functions produced by the Flow with the robot.
# Inside the ipython terminal type:
#
# .. code-block:: ipython
#
# [1] %matplotlib
# [2] robot.plot()
#
# to compare the real and imaginary part of the macroscopic dielectric function for the different displacements.
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_raman_bse.png?raw=true
# :alt: Band structure of Si in the IBZ and along a k-path
#

15
abipy/examples/flows/run_raman_optic.py
View File

@ -40,17 +40,12 @@ def build_flow(options):
if not options.workdir:
options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_","flow_")
manager = options.manager
#shell_manager = manager.to_shell_manager(mpi_procs=1)
#shell_manager = manager.deepcopy()
#ddk_manager = manager.deepcopy()
flow = flowtk.Flow(options.workdir, manager=manager)
flow = flowtk.Flow(options.workdir, manager=options.manager)
# Generate the different shifts to average
ndiv = 1
shift1D = np.arange(1,2*ndiv+1,2)/(2*ndiv)
all_shifts = [[x,y,z] for x in shift1D for y in shift1D for z in shift1D]
shift1D = np.arange(1, 2*ndiv+1, 2) / (2 * ndiv)
all_shifts = [[x, y, z] for x in shift1D for y in shift1D for z in shift1D]
for structure, eta in zip(displaced_structures, etas):
for ishift,shift in enumerate(all_shifts):
@ -108,7 +103,7 @@ def raman_work(structure, pseudos, ngkpt, shiftk):
multi[2 + idir].set_vars(
iscf=-3,
nband=40,
nband=40,
nbdbuf=5,
nstep=1,
nline=0,
@ -158,7 +153,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=False, tight_layout=True)
@flowtk.flow_main

2
abipy/examples/flows/run_relax.py
View File

@ -102,7 +102,7 @@ if os.getenv("GENERATE_SPHINX_GALLERY", False):
__name__ = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()])
build_flow(options).plot_networkx(tight_layout=True)
build_flow(options).plot_networkx(with_edge_labels=True, tight_layout=True)
@flowtk.flow_main

Some files were not shown because too many files have changed in this diff Show More