Add plot_phdispl

This commit is contained in:
Matteo Giantomassi 2018-01-26 17:33:41 +01:00
parent 8fc326c25e
commit 0330ad222f
9 changed files with 254 additions and 37 deletions

View File

@ -16,7 +16,7 @@ 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, add_fig_kwargs, get_ax_fig_plt, get_axarray_fig_plt,
rotate_ticklabels)
rotate_ticklabels, set_visible)
class Robot(NotebookWriter):
@ -949,10 +949,10 @@ Not all entries are sortable (Please select number-like quantities)""" % (self._
for i, (ax, item) in enumerate(zip(ax_list.ravel(), items)):
self.plot_convergence(item, sortby=sortby, hue=hue, ax=ax, fontsize=fontsize,
marker=marker, show=False)
if i != 0 and ax.legend():
ax.legend().set_visible(False)
if i != len(items) - 1 and ax.xaxis.label:
ax.xaxis.label.set_visible(False)
if i != 0:
set_visible(ax, False, "legend")
if i != len(items) - 1:
set_visible(ax, False, "xlabel")
return fig

View File

@ -590,7 +590,10 @@ class Kpoint(SlotPickleMixin):
return np.any(diff < _ATOL_KDIFF)
def __repr__(self):
return "[%+.3f, %+.3f, %+.3f]" % tuple(self.frac_coords)
s = "[%+.3f, %+.3f, %+.3f]" % tuple(self.frac_coords)
if self.name is not None:
s += " %s" % self.name
return s
def __str__(self):
return self.to_string()
@ -598,7 +601,8 @@ class Kpoint(SlotPickleMixin):
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.name is not None:
s += ", name: %s" % self.name
if self._weight is not None: s += ", weight: %.3f" % self.weight
return s

View File

@ -715,7 +715,7 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
return phbands_plotter
def anacompare_phdos(self, nqsmalls, asr=2, chneut=1, dipdip=1, dos_method="tetra", ngqpt=None,
num_cpus=1, stream=sys.stdout):
verbose=0, num_cpus=1, stream=sys.stdout):
"""
Invoke Anaddb to compute Phonon DOS with different q-meshes. The ab-initio dynamical matrix
reported in the DDB_ file will be Fourier-interpolated on the list of q-meshes specified
@ -728,6 +728,7 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
dos_method: Technique for DOS computation in Possible choices: "tetra", "gaussian" or "gaussian:0.001 eV".
In the later case, the value 0.001 eV is used as gaussian broadening
ngqpt: Number of divisions for the ab-initio q-mesh in the DDB file. Auto-detected if None (default)
verbose: Verbosity level.
num_cpus: Number of CPUs (threads) used to parallellize the calculation of the DOSes. Autodetected if None.
stream: File-like object used for printing.
@ -756,7 +757,8 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
else:
# Threads
print("Computing %d phonon DOS with %d threads" % (len(nqsmalls), num_cpus) )
if verbose:
print("Computing %d phonon DOS with %d threads" % (len(nqsmalls), num_cpus) )
phdoses = [None] * len(nqsmalls)
def worker():
@ -790,8 +792,9 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
for i, phdos in enumerate(phdoses[:-1]):
splined_dos = phdos.spline_on_mesh(last_mesh)
abs_diff = (splined_dos - phdoses[-1]).abs()
print(" Delta(Phdos[%d] - Phdos[%d]) / Phdos[%d]: %f" %
(i, len(phdoses)-1, len(phdoses)-1, abs_diff.integral().values[-1]), file=stream)
if verbose:
print(" Delta(Phdos[%d] - Phdos[%d]) / Phdos[%d]: %f" %
(i, len(phdoses)-1, len(phdoses)-1, abs_diff.integral().values[-1]), file=stream)
# Fill the plotter.
plotter = PhononDosPlotter()

View File

@ -23,7 +23,7 @@ from abipy.core.mixins import AbinitNcFile, Has_Structure, Has_PhononBands, Note
from abipy.core.kpoints import Kpoint, Kpath
from abipy.iotools import ETSF_Reader
from abipy.tools import gaussian, duck
from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, set_axlims, get_axarray_fig_plt
from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, set_axlims, get_axarray_fig_plt, set_visible
from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine
from pymatgen.phonon.dos import CompletePhononDos as PmgCompletePhononDos, PhononDos as PmgPhononDos
@ -485,8 +485,7 @@ class PhononBands(object):
def qindex_qpoint(self, qpoint):
"""Returns (qindex, qpoint) from an integer or a qpoint."""
qindex = self.qindex(qpoint)
qpoint = self.qpoints(qindex)
return qindex, qpoint
return qindex, self.qpoints[qindex]
def get_unstable_modes(self, below_mev=-5.0):
"""
@ -1277,7 +1276,7 @@ class PhononBands(object):
axlist: The axes for the bandstructure plot and the DOS plot. If axlist is None, a new figure
is created and the two axes are automatically generated.
width_ratios: Ratio between the width of the bands plots and the DOS plots.
Used if `axlist` is None
Used if ``axlist`` is None
Returns: |matplotlib-Figure|
"""
@ -1321,6 +1320,111 @@ class PhononBands(object):
return fig
@add_fig_kwargs
def plot_phdispl(self, qpoint, cart_dir=None, ax=None, units="eV",
colormap="viridis", hatches="default", fontsize=12, **kwargs):
"""
Plot vertical bars with the contribution of the different atomic types to the phonon displacements
at a given q-point. The contribution is given by the ratio :math:`\dfrac{|d_{type}|^2} {|d|^2}`
where d is the (complex) phonon displacement in cartesian coordinates and d_{type} selects only the
terms associated to the atomic type ``type``.
Args:
qpoint: integer, vector of reduced coordinates or |Kpoint| object.
cart_dir: "x", "y", or "z" to select a particular Cartesian directions.
None if no projection is wanted.
ax: |matplotlib-Axes| or None if a new figure should be created.
units: Units for phonon frequencies. Possible values in ("eV", "meV", "Ha", "cm-1", "Thz").
Case-insensitive.
colormap: Matplotlib colormap used for atom type.
hatches: List of strings with matplotlib hatching patterns. None or empty list to disable hatching.
fontsize: Legend and title fontsize.
Returns: |matplotlib-Figure|
"""
factor = abu.phfactor_ev2units(units)
icart = {"x": 0, "y": 1, "z": 2, None: None}[cart_dir]
iq, qpoint = self.qindex_qpoint(qpoint)
ax, fig, plt = get_ax_fig_plt(ax=ax)
cmap = plt.get_cmap(colormap)
ntypat = self.structure.ntypesp
ax.set_title("qpoint = %s" % repr(qpoint), fontsize=fontsize)
ax.set_xlabel('Frequency %s' % abu.phunit_tag(units))
#ax.set_ylabel(r"Contribution of atomic types to phonon displacement $d_{q,\nu}**2$")
if icart is None:
ax.set_ylabel(r"${|\vec{d}|^2} / {|\vec{d}({type})|^2}$")
else:
ax.set_ylabel(r"${|\vec{d}_{%s}|^2} / {|\vec{d}({type})|^2}$" % cart_dir)
# {Mg: [0], B: [1, 2]}
symbol2indices = {symbol: np.array(self.structure.indices_from_symbol(symbol))
for symbol in self.structure.symbol_set}
#print(symbol2indices)
width, pad = 4, 1
pad = width + pad
xticks, xticklabels = [], []
if hatches == "default":
hatches = ["/", "\\", "'", "|", "-", "+", "x", "o", "O", ".", "*"]
else:
hatches = list_strings(hatches) if hatches is not None else []
x = 0
for nu in self.branches:
w_qnu = self.phfreqs[iq, nu] * factor
dcart_qnu = np.reshape(self.phdispl_cart[iq, nu], (len(self.structure), 3))
d2norm = sum(np.linalg.norm(d)**2 for d in dcart_qnu)
# Make a bar plot with rectangles bounded by (x - width/2, x + width/2, bottom, bottom + height)
# The align keyword controls if x is interpreted as the center or the left edge of the rectangle.
bottom, height = 0.0, 0.0
for itype, (symbol, inds) in enumerate(symbol2indices.items()):
if icart is None:
assert all(len(d) == 3 for d in dcart_qnu[inds])
height = sum(np.linalg.norm(d)**2 for d in dcart_qnu[inds]) / d2norm
else:
height = sum(np.linalg.norm(d)**2 for d in dcart_qnu[inds, icart]) / d2norm
ax.bar(x, height, width, bottom, align="center",
color=cmap(float(itype) / max(1, ntypat - 1)),
label=symbol if nu == 0 else None, edgecolor='black',
hatch=hatches[itype % len(hatches)] if hatches else None,
)
bottom += height
xticks.append(x)
xticklabels.append("%.3f" % w_qnu)
x += (width + pad) / 2
ax.set_xticks(xticks)
ax.set_xticklabels((xticklabels))
ax.legend(loc="best", fontsize=fontsize, shadow=False)
#ax.legend(bbox_to_anchor=(0, 1.02, 1, 0.2), loc="lower left", mode="expand", borderaxespad=0, ncol=ntypat)
return fig
@add_fig_kwargs
def plot_phdispl_cartdirs(self, qpoint, units="eV", colormap="viridis", hatches="default", fontsize=12, **kwargs):
"""
Plot vertical bars with the contribution of the different atomic types to the phonon displacements
at a gven q-point. See plot_phdispl for API documentation.
"""
ax_list, fig, plt = get_axarray_fig_plt(None, nrows=3, ncols=1,
sharex=True, sharey=True, squeeze=False)
cart_dirs = ["x", "y", "z"]
for i, (cart_dir, ax) in enumerate(zip(cart_dirs, ax_list.ravel())):
self.plot_phdispl(qpoint, cart_dir=cart_dir, ax=ax, units=units, colormap=colormap,
fontsize=fontsize, hatches=hatches, show=False)
# Disable artists.
if i != 0:
set_visible(ax, False, "legend", "title")
if i != len(cart_dirs) - 1:
set_visible(ax, False, "xlabel")
return fig
def get_dataframe(self):
"""
Return a |pandas-DataFrame| with the following columns:
@ -1478,7 +1582,7 @@ class PhononBands(object):
the maximum breaking with sign
the absolute value of the maximum breaking
"""
gamma_ind = self.qpoints.index((0,0,0))
gamma_ind = self.qpoints.index((0, 0, 0))
ind = self.acoustic_indices(gamma_ind, threshold=threshold, raise_on_no_indices=raise_on_no_indices)
asr_break = self.phfreqs[0, ind] * abu.phfactor_ev2units(units)
@ -2071,6 +2175,7 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
"""
Human-readable string with useful information such as structure...
Args:
verbose: Verbosity level.
"""
lines = []; app = lines.append
@ -2078,7 +2183,6 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
app(marquee("File Info", mark="="))
app(self.filestat(as_string=True))
app("")
app(self.structure.to_string(verbose=verbose, title="Structure"))
app("")
@ -2867,7 +2971,7 @@ class PhononBandsPlotter(NotebookWriter):
mode_range: Only bands such as ``mode_range[0] <= nu_index < mode_range[1]`` are included in the plot.
units: Units for phonon plots. Possible values in ("eV", "meV", "Ha", "cm-1", "Thz"). Case-insensitive.
swarm: True to show the datapoints on top of the boxes
ax: matplotlib :class:`Axes` or None if a new figure should be created.
ax: |matplotlib-Axes| or None if a new figure should be created.
kwargs: Keyword arguments passed to seaborn_ boxplot.
"""
frames = []

View File

@ -40,6 +40,7 @@ class PhononBandsTest(AbipyTest):
assert any(q.name is not None for q in nc.qpoints)
same_phbands_nc = PhononBands.as_phbands(nc)
self.assert_equal(same_phbands_nc.phfreqs, phbands.phfreqs)
assert phbands.phdispl_cart.shape == (phbands.nqpt, phbands.num_branches, phbands.num_branches)
# a + b gives plotter
assert hasattr(same_phbands_nc + phbands, "combiplot")
@ -103,6 +104,8 @@ class PhononBandsTest(AbipyTest):
assert phbands.plot_fatbands(units="ha", qlabels=qlabels, show=False)
assert phbands.plot_fatbands(phdos_file=abidata.ref_file("trf2_5.out_PHDOS.nc"), units="thz", show=False)
assert phbands.plot_colored_matched(units="cm^-1", show=False)
assert phbands.plot_phdispl(qpoint=(0, 0, 0), units="cm^-1", hatches=None, show=False)
assert phbands.plot_phdispl_cartdirs(qpoint=(0, 0, 0), units="cm^-1", show=False)
assert phbands.boxplot(units="ev", mode_range=[2, 4], show=False)
# Cannot compute PHDOS with q-path
@ -153,8 +156,9 @@ class PhbstFileTest(AbipyTest):
ii = ncfile.qindex(qpt)
#print("iq", iq, "qpt", qpt, "ii", ii, "qpoints[ii]", ncfile.qpoints[ii])
same_ii, same_qpt = ncfile.qindex_qpoint(ii)
assert same_ii == ii
assert qpt == same_qpt
assert same_ii == ii and qpt == same_qpt
same_ii, same_qpt = ncfile.phbands.qindex_qpoint(qpt)
assert same_ii == ii and qpt == same_qpt
qpoint = ncfile.qpoints[0]
frame = ncfile.get_phframe(qpoint)

View File

@ -184,6 +184,10 @@ class GsrFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, Notebo
"""
return self.reader.read_abinit_xcfunc()
@lazy_property
def energy_terms(self):
return self.reader.read_energy_terms(unit="eV")
@lazy_property
def params(self):
""":class:`OrderedDict` with parameters that might be subject to convergence studies."""
@ -273,7 +277,9 @@ if gsr.ebands.kpoints.is_ibz:
class EnergyTerms(AttrDict):
"""Contributions to the total GS energy. See energies_type in m_energies.F90"""
"""
Contributions to the total GS energy. See energies_type in m_energies.F90.
"""
_NAME2DOC = OrderedDict([
# (Name, help)
@ -317,6 +323,15 @@ class EnergyTerms(AttrDict):
return self.to_string(with_doc=False)
__repr__ = __str__
def to_string(self, verbose=0, with_doc=True):
"""String representation, with documentation if with_doc."""
lines = [str(self.table)]
if with_doc:
for k, doc in self._NAME2DOC.items():
lines.append("%s: %s" % (k, doc))
return "\n".join(lines)
@property
def table(self):
"""PrettyTable object with the results."""
@ -326,14 +341,13 @@ class EnergyTerms(AttrDict):
table.add_row([k, self[k]])
return table
def to_string(self, verbose=0, with_doc=True):
"""String representation, with documentation if with_doc."""
lines = [str(self.table)]
if with_doc:
for k, doc in self._NAME2DOC.items():
lines.append("%s: %s" % (k, doc))
return "\n".join(lines)
#def get_dataframe(self):
# """
# Return a |pandas-DataFrame|
# """
# d = {k: float(self[k]) for k in self}
# return pd.DataFrame(d, index=[None], columns=list(d.keys()))
# #return pd.DataFrame(d, columns=list(d.keys()))
class GsrReader(ElectronsReader):
@ -486,6 +500,69 @@ class GsrRobot(Robot, RobotWithEbands):
dataframe = pd.DataFrame(rows, index=index, columns=list(rows[0].keys()) if rows else None)
return dict2namedtuple(fits=fits, dataframe=dataframe)
def get_energyterms_dataframe(self, abspath=False):
rows, row_names = [], []
for label, gsr in self.items():
row_names.append(label)
rows.append(gsr.energy_terms)
row_names = row_names if not abspath else self._to_relpaths(row_names)
df = pd.DataFrame(rows, index=row_names, columns=list(rows[0].keys()))
first_row = df.iloc[[0]].values[0]
df = df.apply(lambda row: row - first_row, axis=1)
return df
@add_fig_kwargs
def plot_energyterms(self, fontsize=6, **kwargs):
items = ["e_eigenvalues", "e_localpsp", "e_ewald"]
# Build grid of plots.
nrows, ncols = len(items), 1
ax_list, fig, plt = get_axarray_fig_plt(None, nrows=nrows, ncols=ncols,
sharex=True, sharey=False, squeeze=False)
df = self.get_energyterms_dataframe()
#for i, (item, ax) in enumerate(zip(items, ax_list.ravel())):
# print(item, "\n", df[item])
# df.plot(x=None, y=item, kind='bar', ax=ax, legend=i==0, fontsize=fontsize)
for i, (item, ax) in enumerate(zip(items, ax_list.ravel())):
#values = [float(afile.energy_terms[item]) for afile in self.abifiles]
bottom = 0
width, pad = 4, 1
pad = width + pad
xticks, xticklabels = [], []
x = 0
for j, (label, abifile) in enumerate(self):
height = float(abifile.energy_terms[item])
if j == 0: h0 = height
height -= h0
ax.bar(x, height, width, bottom, align="center",
label=label,
#color=cmap(float(itype) / max(1, ntypat - 1)),
edgecolor='black',
#hatch=hatches[itype % len(hatches)] if hatches else None,
)
xticks.append(x)
xticklabels.append(label)
x += (width + pad) / 2
ax.set_xticks(xticks)
ax.set_xticklabels((xticklabels))
ax.legend(loc="best", fontsize=fontsize, shadow=False)
#rows, row_names = [], []
#for label, gsr in self.items():
# row_names.append(label)
# rows.append(gsr.energy_terms)
#row_names = row_names if not abspath else self._to_relpaths(row_names)
#data = pd.DataFrame(rows, index=row_names, columns=list(rows[0].keys()))
#ax, fig, plt = get_ax_fig_plt(ax=None)
#sns.barplot(x="day", y="total_bill", hue="sex", data=data)
return fig
@add_fig_kwargs
def gridplot_eos(self, eos_names="all", fontsize=6, **kwargs):
"""

View File

@ -107,18 +107,17 @@ if __name__ == "__main__":
#
# abicomp.py ddb flow_mgb2_phonons_nkpt_tsmear/w*/outdata/*_DDB -ipy
#
# to build a robot from all the output DDB files and start the ipython shell.
# to build a robot from the output DDB files and start the ipython shell.
#
# then, inside the ipython shell, use:
#
# .. code-block:: ipython
#
# In [1]: %matplotlib
# In [2]: p = robot.anaget_phonon_plotters()
# In [3]: p.phbands_plotter.combiplot()
# In [2]: r = robot.anaget_phonon_plotters(nqsmall=0)
# In [3]: r.phbands_plotter.gridplot_with_hue("tsmear")
#
# to compute the phonon bands and the DOS for all DDB files with Anaddb
# and plot the results on the same figure.
# to compute the phonon bands with Anaddb and plot the results groupe by "tsmear".
#
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_mgb2_edoses.png?raw=true
# .. image:: https://github.com/abinit/abipy_assets/blob/master/run_mgb2_phonons_nkpt_tsmear.png?raw=true
# :alt: Convergence of electronic DOS in MgB2 wrt k-points.

View File

@ -728,8 +728,12 @@ Usage example:
TIP:
Use Unix find to select all files with the a given extension and pass them to abicomp.py.
For instance:
The python code operates on a list of files/directories passed via the command line interface.
The arguments are interpreted by the shell before invoking the script.
This means that one can use the bash syntax and the unix tools to precompute the list of files/directories.
For example, one can use Unix find to select all files with the a given extension and pass them to abicomp.py.
For command:
abicomp.py structure `find . -name "*_GSR.nc"`
@ -744,6 +748,14 @@ and this trick can be used to select files/directories as in:
abicomp.py structure w1/t{2..4}/outdata/*_GSR.nc
The [!name] syntax can be used to exclude patterns, so
abicomp.py structure w1/t[!2]*/outdata/*_GSR.nc
excludes all the GSR.nc files in the t2/outdata directory.
See also http://wiki.bash-hackers.org/syntax/pattern
NOTE: The `gsr`, `ddb`, `sigres`, `mdf` commands use robots to analyze files.
In this case, one can provide a list of files and/or list of directories on the command-line interface e.g.:

View File

@ -60,6 +60,20 @@ def set_axlims(ax, lims, axname):
return left, right
def set_visible(ax, boolean, *args):
"""
Hide/Show the artists of axis ax listed in args.
"""
if "legend" in args and ax.legend():
ax.legend().set_visible(boolean)
if "title" in args and ax.title:
ax.title.set_visible(boolean)
if "xlabel" in args and ax.xaxis.label:
ax.xaxis.label.set_visible(boolean)
if "ylabel" in args and ax.yaxis.label:
ax.yaxis.label.set_visible(boolean)
def rotate_ticklabels(ax, rotation, axname="x"):
"""Rotate the ticklables of axis ``ax``"""
if "x" in axname: