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

Merge branch 'qha' of https://github.com/gpetretto/abipy into gp_qha 

* 'qha' of https://github.com/gpetretto/abipy:
  use amu_symbol
  finite difference gruneisen
This commit is contained in:
Matteo Giantomassi 2018-08-09 21:23:47 +02:00
parent c4b4854f48 ccc46b1bdd
commit 297ae2d758
3 changed files with 215 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

211
abipy/dfpt/gruneisen.py
View File

@ -13,16 +13,22 @@ from monty.termcolor import cprint
from monty.collections import AttrDict
from monty.functools import lazy_property
from abipy.core.kpoints import Kpath, IrredZone, KSamplingInfo
from abipy.core.mixins import AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter
from abipy.core.mixins import AbinitNcFile, Has_Structure, NotebookWriter
from abipy.abio.inputs import AnaddbInput
from abipy.dfpt.phonons import PhononBands, PhononBandsPlotter, PhononDos
from abipy.dfpt.phonons import PhononBands, PhononBandsPlotter, PhononDos, match_eigenvectors, get_dyn_mat_eigenvec
from abipy.dfpt.ddb import DdbFile
from abipy.iotools import ETSF_Reader
from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, get_axarray_fig_plt, set_axlims
from abipy.flowtk import AnaddbTask
from abipy.tools.derivatives import finite_diff
#from abipy.tools import duck
from abipy.core.func1d import Function1D
from pymatgen.core.units import amu_to_kg, bohr_to_ang
from pymatgen.core.units import amu_to_kg
from pymatgen.core.periodic_table import Element
try:
from functools import lru_cache
except ImportError: # py2k
from abipy.tools.functools_lru_cache import lru_cache
# DOS name --> meta-data
@ -109,12 +115,16 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
@lazy_property
def wvols_qibz(self):
"""Phonon frequencies on reagular grid for the different volumes in eV """
return self.reader.read_value("gruns_wvols_qibz") * abu.Ha_eV
w = self.reader.read_value("gruns_wvols_qibz", default=None)
if w is None:
return None
else:
return w * abu.Ha_eV
@lazy_property
def qibz(self):
"""q-points in the irreducible brillouin zone"""
return self.reader.read_value("gruns_qibz")
return self.reader.read_value("gruns_qibz", default=None)
@lazy_property
def gvals_qibz(self):
@ -130,6 +140,34 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
"""List of |PhononBands| objects corresponding to the different volumes."""
return self.reader.read_phbands_on_qpath()
@lazy_property
def structures(self):
"""List of structures"""
return self.reader.read_structures()
@lazy_property
def volumes(self):
"""List of volumes"""
return [s.volume for s in self.structures]
@lazy_property
def phdispl_cart_qibz(self):
"""Eigendisplacements for the modes on the qibz"""
return self.reader.read_value("gruns_phdispl_cart_qibz", cmode="c")
@property
def nvols(self):
"""Number of volumes"""
return len(self.structures)
@lazy_property
def amu_symbol(self):
"""Atomic mass units"""
amu_list = self.reader.read_value("atomic_mass_units")
atomic_numbers = self.reader.read_value("atomic_numbers")
amu = {Element.from_Z(at).symbol: a for at, a in zip(atomic_numbers, amu_list)}
return amu
def to_dataframe(self):
"""
Return a |pandas-DataFrame| with the following columns:
@ -233,8 +271,8 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
showing the value and the sign of the Grunesein parameters.
Args:
fill_with: Define the quantity used to plot stripes. "gruns" for Grunesein parameters,
"groupv" for phonon group velocities.
fill_with: Define the quantity used to plot stripes. "gruns" for Grunesein parameters, "gruns_fd" for
Grunesein parameters calculated with finite differences, "groupv" for phonon group velocities.
gamma_fact: Scaling factor for Grunesein parameters.
Up triangle for positive values, down triangles for negative values.
alpha: The alpha blending value for the markers between 0 (transparent) and 1 (opaque).
@ -277,11 +315,15 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
if fill_with == "gruns":
max_gamma = np.abs(phbands.grun_vals).max()
values = phbands.grun_vals
elif fill_with == "groupv":
# TODO: units?
dwdq_qpath = self.reader.read_value("gruns_dwdq_qpath")
groupv = np.linalg.norm(dwdq_qpath, axis=-1)
max_gamma = np.abs(groupv).max()
values = np.linalg.norm(dwdq_qpath, axis=-1)
max_gamma = np.abs(values).max()
elif fill_with == "gruns_fd":
max_gamma = np.abs(self.grun_vals_finite_differences(match_eigv=True)).max()
values = self.grun_vals_finite_differences(match_eigv=True)
else:
raise ValueError("Unsupported fill_with: `%s`" % fill_with)
@ -298,16 +340,16 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
for nu in branch_range:
omegas = phbands.phfreqs[:, nu].copy() * factor
if fill_with == "gruns":
if fill_with.startswith("gruns"):
# Must handle positive-negative values
sizes = phbands.grun_vals[:, nu].copy() * (gamma_fact * 0.02 * max_omega / max_gamma)
sizes = values[:, nu].copy() * (gamma_fact * 0.02 * max_omega / max_gamma)
yup = omegas + np.where(sizes >= 0, sizes, 0)
ydown = omegas + np.where(sizes < 0, sizes, 0)
ax_bands.fill_between(xvals, omegas, yup, alpha=alpha, facecolor="red")
ax_bands.fill_between(xvals, ydown, omegas, alpha=alpha, facecolor="blue")
elif fill_with == "groupv":
sizes = groupv[:, nu].copy() * (gamma_fact * 0.04 * max_omega / max_gamma)
sizes = values[:, nu].copy() * (gamma_fact * 0.04 * max_omega / max_gamma)
ydown, yup = omegas - sizes / 2, omegas + sizes / 2
ax_bands.fill_between(xvals, ydown, yup, alpha=alpha, facecolor="red")
@ -348,8 +390,8 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
of the phonon frequencies.
Args:
values: Define the plotted quantity. "gruns" for Grunesein parameters,
"groupv" for phonon group velocities.
values: Define the plotted quantity. "gruns" for Grunesein parameters, "gruns_fd" for Grunesein
parameters calculated with finite differences, "groupv" for phonon group velocities.
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.
@ -363,6 +405,8 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
elif values == "groupv":
# TODO: units?
y = np.linalg.norm(self.reader.read_value("gruns_dwdq_qibz"), axis=-1)
elif values == "gruns_fd":
y = self.gvals_qibz_finite_differences(match_eigv=True)
else:
raise ValueError("Unsupported values: `%s`" % values)
@ -375,7 +419,7 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
ax.scatter(w.flatten(), y.flatten(), **kwargs)
ax.set_xlabel('Frequency %s' % abu.phunit_tag(units))
if values == "gruns":
if values.startswith("gruns"):
ax.set_ylabel('Gruneisen')
elif values == "groupv":
ax.set_ylabel('|v|')
@ -425,8 +469,8 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
and match_bands=True to obtained reasonable paths.
Args:
values: Define the plotted quantity. "gruns" for Grunesein parameters,
"groupv" for phonon group velocities.
values: Define the plotted quantity. "gruns" for Grunesein parameters, "gruns_fd" for Grunesein
parameters calculated with finite differences, "groupv" for phonon group velocities.
ax: |matplotlib-Axes| or None if a new figure should be created.
branch_range: Tuple specifying the minimum and maximum branch index to plot (default: all
branches are plotted).
@ -444,6 +488,8 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
elif values == "groupv":
# TODO: units?
y = np.linalg.norm(self.split_dwdq, axis=-1)
elif values == "gruns_fd":
y = self.split_gruns_finite_differences(match_eigv=True)
else:
raise ValueError("Unsupported values: `%s`" % values)
@ -458,7 +504,7 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
ax, fig, plt = get_ax_fig_plt(ax=ax)
phbands.decorate_ax(ax, units=None, qlabels=qlabels)
if values == "gruns":
if values in ("gruns", "gruns_fd"):
ax.set_ylabel('Gruneisen')
elif values == "groupv":
ax.set_ylabel('|v|')
@ -708,6 +754,64 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
return gruns
@lru_cache()
def grun_vals_finite_differences(self, match_eigv=True):
"""
Gruneisen parameters along the high symmetry path calculated with finite differences.
"""
if not self.phbands_qpath_vol:
raise ValueError("Finite differences require the phonon bands")
phbands = np.array([b.phfreqs for b in self.phbands_qpath_vol])
if match_eigv:
eig = np.array([b.dyn_mat_eigenvect for b in self.phbands_qpath_vol])
else:
eig = None
dv = np.abs(self.volumes[0] - self.volumes[1])
return calculate_gruns_finite_differences(phbands, eig, self.iv0, self.structure.volume, dv)
@lru_cache()
def split_gruns_finite_differences(self, match_eigv=True):
"""
Splits the values of the finite differences gruneisen along a path like for the phonon bands
"""
try:
return self._split_gruns_fd
except AttributeError:
# trigger the generation of the split in the phbands
self.phbands_qpath_vol[self.iv0].split_phfreqs
indices = self.phbands_qpath_vol[self.iv0]._split_indices
g = self.grun_vals_finite_differences(match_eigv=match_eigv)
self._split_gruns_fd = [np.array(g[indices[i]:indices[i + 1] + 1]) for i in range(len(indices) - 1)]
return self._split_gruns_fd
@lru_cache()
def gvals_qibz_finite_differences(self, match_eigv=True):
"""
Gruneisen parameters in the irreducible brillouin zone calculated with finite differences.
"""
if self.wvols_qibz is None:
raise ValueError("Finite differences require wvols_qibz")
if match_eigv:
eig = np.zeros_like(self.phdispl_cart_qibz)
for i in range(self.nvols):
eig[:, i] = get_dyn_mat_eigenvec(self.phdispl_cart_qibz[:, i], self.structures[i],
amu_symbol=self.amu_symbol)
eig = eig.transpose((1, 0, 2, 3))
else:
eig = None
dv = np.abs(self.volumes[0] - self.volumes[1])
return calculate_gruns_finite_differences(self.wvols_qibz.transpose(1, 0, 2), eig, self.iv0,
self.structure.volume, dv)
class GrunsReader(ETSF_Reader):
"""
@ -806,18 +910,12 @@ class GrunsReader(ETSF_Reader):
lattices = self.read_value("gruns_rprimd") * abu.Bohr_Ang #, "dp", "three, three, gruns_nvols")
gruns_xred = self.read_value("gruns_xred") #, "dp", "three, number_of_atoms, gruns_nvols")
structures = self.read_structures()
phbands_qpath_vol = []
for ivol in range(self.num_volumes):
# TODO structure depends on:
# volumes
# non_anal_ph
# amu: DONE
# phdispl_cart DONE
if ivol == self.iv0:
structure = self.structure
else:
structure = self.structure.__class__(lattices[ivol], self.structure.species, gruns_xred[ivol])
structure = structures[ivol]
# TODO non_anal_ph
qpoints = Kpath(structure.reciprocal_lattice, qfrac_coords)
phdispl_cart = phdispl_cart_qptsvol[:, ivol].copy()
phb = PhononBands(structure, qpoints, freqs_vol[:, ivol], phdispl_cart, non_anal_ph=None, amu=amuz)
@ -841,3 +939,58 @@ class GrunsReader(ETSF_Reader):
amuz[znucl_typat[type_idx]] = amu_typat[type_idx]
return amuz
def read_structures(self):
"""
Resturns a list of structures at the different volumes
"""
lattices = self.read_value("gruns_rprimd") * abu.Bohr_Ang # , "dp", "three, three, gruns_nvols")
gruns_xred = self.read_value("gruns_xred") # , "dp", "three, number_of_atoms, gruns_nvols")
structures = []
for ivol in range(self.num_volumes):
if ivol == self.iv0:
structures.append(self.structure)
else:
structures.append(self.structure.__class__(lattices[ivol], self.structure.species, gruns_xred[ivol]))
return structures
def calculate_gruns_finite_differences(phfreqs, eig, iv0, volume, dv):
"""
Calculates the Gruneisen parameters from finite differences on the phonon frequencies. Uses the eigenvectors
to match the frequencies at different volumes.
Args:
phfreqs: numpy array with the phonon frequencies at different volumes. Shape (nvols, nqpts, 3*natoms)
eig: numpy array with the phonon eigenvectors at the different volumes. Shape (nvols, nqpts, 3*natoms, 3*natoms)
If None simple ordering will be used.
iv0: index of the 0 volume.
volume: volume of the structure at the central volume.
dv: volume variation.
Returns:
A numpy array with the gruneisen parameters. Shape (nqpts, 3*natoms)
"""
phfreqs = phfreqs.copy()
nvols = phfreqs.shape[0]
if eig is not None:
for i in range(nvols):
if i == iv0:
continue
for iq in range(phfreqs.shape[1]):
ind = match_eigenvectors(eig[iv0, iq], eig[i, iq])
phfreqs[i, iq] = phfreqs[i, iq][ind]
acc = nvols - 1
g = np.zeros_like(phfreqs[0])
for iq in range(phfreqs.shape[1]):
for im in range(phfreqs.shape[2]):
w = phfreqs[iv0, iq, im]
if w == 0:
g[iq, im] = 0
else:
g[iq, im] = - finite_diff(phfreqs[:, iq, im], dv, order=1, acc=acc)[iv0] * volume / w
return g

21
abipy/dfpt/phonons.py
View File

@ -4057,7 +4057,7 @@ class InteratomicForceConstants(Has_Structure):
# TODO: amu should become mandatory.
def get_dyn_mat_eigenvec(phdispl, structure, amu=None):
def get_dyn_mat_eigenvec(phdispl, structure, amu=None, amu_symbol=None):
"""
Converts the phonon displacements to the orthonormal eigenvectors of the dynamical matrix.
Small discrepancies with the original values may be expected due to the different values of the atomic masses in
@ -4075,20 +4075,29 @@ def get_dyn_mat_eigenvec(phdispl, structure, amu=None):
should match the q points.
structure: |Structure| object.
amu: dictionary that associates the atomic numbers present in the structure to the values of the atomic
mass units used for the calculation. If None, values from pymatgen will be used. Note that this will
almost always lead to inaccuracies in the conversion.
mass units used for the calculation. Incompatible with amu_sumbol. If None and amu_symbol is None, values
from pymatgen will be used. Note that this will almost always lead to inaccuracies in the conversion.
amu_symbol: dictionary that associates the symbol present in the structure to the values of the atomic
mass units used for the calculation. Incompatible with amu. If None and amu_symbol is None, values from
pymatgen will be used. that this will almost always lead to inaccuracies in the conversion.
Returns:
A |numpy-array| of the same shape as phdispl containing the eigenvectors of the dynamical matrix
"""
eigvec = np.zeros(np.shape(phdispl), dtype=np.complex)
if amu is None:
if amu is not None and amu_symbol is not None:
raise ValueError("Only one between amu and amu_symbol should be provided!")
if amu is not None:
amu_symbol = {Element.from_Z(n).symbol: v for n, v in amu.items()}
if amu_symbol is None:
warnings.warn("get_dyn_mat_eigenvec has been called with amu=None. Eigenvectors may not be orthonormal.")
amu = {e.number: e.atomic_mass for e in structure.composition.elements}
amu_symbol = {e.symbol: e.atomic_mass for e in structure.composition.elements}
for j, a in enumerate(structure):
eigvec[...,3*j:3*(j+1)] = phdispl[...,3*j:3*(j+1)] * np.sqrt(amu[a.specie.number]*abu.amu_emass) / abu.Bohr_Ang
eigvec[...,3*j:3*(j+1)] = phdispl[...,3*j:3*(j+1)] * np.sqrt(amu_symbol[a.specie.symbol]*abu.amu_emass) / abu.Bohr_Ang
return eigvec

19
abipy/dfpt/tests/test_gruneisen.py
View File

@ -7,7 +7,7 @@ import abipy.data as abidata
from abipy.core.testing import AbipyTest
from abipy import abilab
from abipy.dfpt.gruneisen import GrunsNcFile
from abipy.dfpt.gruneisen import GrunsNcFile, calculate_gruns_finite_differences
class GrunsFileTest(AbipyTest):
@ -42,6 +42,9 @@ class GrunsFileTest(AbipyTest):
self.assertAlmostEqual(ncfile.debye_temp, 429.05702577371898)
self.assertAlmostEqual(ncfile.acoustic_debye_temp, 297.49152615955893)
ncfile.grun_vals_finite_differences(match_eigv=False)
ncfile.gvals_qibz_finite_differences(match_eigv=False)
if self.has_matplotlib():
assert ncfile.plot_doses(title="DOSes", show=False)
assert ncfile.plot_doses(with_idos=False, xlims=None, show=False)
@ -50,6 +53,7 @@ class GrunsFileTest(AbipyTest):
assert ncfile.plot_phbands_with_gruns(title="bands with gamma markers + DOSes", show=False)
assert ncfile.plot_phbands_with_gruns(with_doses=None, gamma_fact=2, units="cm-1", show=False)
assert ncfile.plot_phbands_with_gruns(fill_with="groupv", gamma_fact=2, units="cm-1", show=False)
assert ncfile.plot_phbands_with_gruns(fill_with="gruns_fd", gamma_fact=2, units="cm-1", show=False)
plotter = ncfile.get_plotter()
assert plotter.combiboxplot(show=False)
@ -57,9 +61,11 @@ class GrunsFileTest(AbipyTest):
assert ncfile.plot_gruns_scatter(units='cm-1', show=False)
assert ncfile.plot_gruns_scatter(values="groupv", units='cm-1', show=False)
assert ncfile.plot_gruns_scatter(values="gruns_fd", units='cm-1', show=False)
assert ncfile.plot_gruns_bs(match_bands=True, show=False)
assert ncfile.plot_gruns_bs(values="groupv", match_bands=False, show=False)
assert ncfile.plot_gruns_bs(values="gruns_fd", match_bands=False, show=False)
if self.has_nbformat():
assert ncfile.write_notebook(nbpath=self.get_tmpname(text=True))
@ -73,3 +79,14 @@ class GrunsFileTest(AbipyTest):
ddb_list = [os.path.join(path, "mp-149_{:+d}_DDB".format(s)) for s in strains]
g = GrunsNcFile.from_ddb_list(ddb_list, ndivsm=3, nqsmall=3)
class FunctionsTest(AbipyTest):
def test_calculate_gruns_finite_differences(self):
phfreqs = np.array([[[0, 0, 0]], [[1, 2, 3]], [[2, 6, 4]]])
eig = np.array([[[[1, 0, 0], [0, 1, 0], [0, 0, 1]]], [[[1, 0, 0], [0, 0, 1], [0, 1, 0]]],
[[[1, 0, 0], [0, 1, 0], [0, 0, 1]]]])
g = calculate_gruns_finite_differences(phfreqs, eig, iv0=1, volume=1, dv=1)
self.assertArrayEqual(g, [[-1, -1, -1]])