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

Use pmg Tensor for anaddbnc.emacro

This commit is contained in:
Matteo Giantomassi 2018-08-08 17:48:49 +02:00
parent a4f35e1302
commit 61b6133efd
5 changed files with 194 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
abipy/core/structure.py
View File

@ -1001,7 +1001,6 @@ class Structure(pymatgen.Structure, NotebookWriter):
return None
def get_symbol2indices(self):
"""
Return a dictionary mapping chemical symbols to numpy array with the position of the atoms.

3
abipy/core/testing.py
View File

@ -570,7 +570,6 @@ def notebook_run(path):
nb = nbformat.read(fout, nbformat.current_nbformat)
errors = [output for cell in nb.cells if "outputs" in cell
for output in cell["outputs"]\
if output.output_type == "error"]
for output in cell["outputs"] if output.output_type == "error"]
return nb, errors

253
abipy/dfpt/ddb.py
View File

@ -16,6 +16,8 @@ from monty.string import marquee, list_strings
from monty.collections import AttrDict, dict2namedtuple, tree
from monty.functools import lazy_property
from monty.termcolor import cprint
from pymatgen.analysis.elasticity import Tensor, Stress, ElasticTensor
from pymatgen.core.units import eV_to_Ha, bohr_to_angstrom, Energy
from abipy.flowtk import NetcdfReader, AnaddbTask
from abipy.core.mixins import TextFile, Has_Structure, NotebookWriter
from abipy.core.symmetries import AbinitSpaceGroup
@ -27,15 +29,10 @@ from abipy.dfpt.phonons import PhononDosPlotter, PhononBandsPlotter, Interatomic
from abipy.dfpt.tensors import DielectricTensor
from abipy.dfpt.elastic import ElasticData
from abipy.core.abinit_units import phfactor_ev2units, phunit_tag #Ha_cmm1,
from pymatgen.analysis.elasticity import Stress, ElasticTensor
from pymatgen.core.units import eV_to_Ha, bohr_to_angstrom, Energy
from abipy.tools.plotting import Marker, add_fig_kwargs, get_ax_fig_plt, set_axlims
from abipy.tools import duck
from abipy.abio.robots import Robot
import logging
logger = logging.getLogger(__name__)
try:
from functools import lru_cache
except ImportError: # py2k
@ -1513,7 +1510,7 @@ class Becs(Has_Structure):
def __init__(self, becs_arr, structure, chneut, order="c"):
"""
Args:
becs_arr: (3, 3, natom) array with the Born effective charges in Cartesian coordinates.
becs_arr: [3, 3, natom] array with the Born effective charges in Cartesian coordinates.
structure: |Structure| object.
chneut: Option used for the treatment of the Charge Neutrality requirement
for the effective charges (anaddb input variable)
@ -1523,12 +1520,15 @@ class Becs(Has_Structure):
self._structure = structure
self.chneut = chneut
# Values is a numpy array while zstars is a list of Tensor objects.
self.values = np.empty((len(structure), 3, 3))
for i, bec in enumerate(becs_arr):
mat = becs_arr[i]
if order.lower() == "f": mat = mat.T.copy()
self.values[i] = mat
self.zstars = [Tensor(mat) for mat in self.values]
@property
def structure(self):
"""|Structure| object."""
@ -1539,13 +1539,17 @@ class Becs(Has_Structure):
def to_string(self, verbose=0):
"""String representation."""
lines = []
app = lines.append
app("Born effective charges computed with chneut: %d\n" % self.chneut)
for site, bec in zip(self.structure, self.values):
app("Z* at site: %s" % repr(site))
app(str(bec))
app("")
lines = []; app = lines.append
app("Born effective charges in Cartesian coordinates (Voigt notation)")
app(self.get_voigt_dataframe().to_string())
app("")
if verbose:
app("Born effective charges (full tensor)")
for site, bec in zip(self.structure, self.values):
app("Z* at site: %s" % repr(site))
app(str(bec))
app("")
# Add info on the bec sum rule.
app("Born effective charge neutrality sum-rule with chneut: %d\n" % self.chneut)
@ -1555,13 +1559,35 @@ class Becs(Has_Structure):
@property
def sumrule(self):
"""Born effective charge neutrality sum-rule."""
"""[3, 3] matrix with Born effective charge neutrality sum-rule."""
return self.values.sum(axis=0)
# TODO: Deprecated
def check_sumrule(self, stream=sys.stdout):
stream.write("Born effective charge neutrality sum-rule with chneut: %d\n" % self.chneut)
stream.write(str(self.sumrule))
def _repr_html_(self):
"""Integration with jupyter notebooks."""
return self.get_voigt_dataframe()._repr_html_()
def get_voigt_dataframe(self, select_symbols=None):
"""
Return |pandas-DataFrame| with Voigt indices as columns and natom rows.
Args:
select_symbols: String or list of strings with chemical symbols.
Used to select only atoms of this type.
"""
select_symbols = set(list_strings(select_symbols)) if select_symbols is not None else None
columns = ["xx", "yy", "zz", "yz", "xz", "xy"]
rows = []
for site, zstar in zip(self.structure, self.zstars):
if select_symbols is not None and site.specie.symbol not in select_symbols: continue
d = OrderedDict()
d["element"] = site.specie.symbol
d["frac_coords"] = site.frac_coords
for k, v in zip(columns, zstar.voigt):
d[k] = v
rows.append(d)
return pd.DataFrame(rows, index=None, columns=list(rows[0].keys()))
class DielectricTensorGenerator(Has_Structure):
@ -1636,7 +1662,7 @@ class DielectricTensorGenerator(Has_Structure):
phfreqs = phbands.phfreqs[gamma_index]
emacro = anaddbnc.emacro.cartesian_tensor
emacro = anaddbnc.emacro
oscillator_strength = anaddbnc.oscillator_strength
return cls(phfreqs, oscillator_strength, emacro, anaddbnc.structure)
@ -1785,73 +1811,70 @@ class DdbRobot(Robot):
#
# return np.array(qpoints)
# Uncomment this method to use `abicomp ddb */out_DDB -p`
#def to_string(self, **kwargs):
# self.find_duplicated_entries()
# DEBUGGING CODE (do not remove)
#def find_duplicated_entries(self, std_tol=1e-5, verbose=1):
# """
# Check for duplicated entries in the list of ddb files
def find_duplicated_entries(self, std_tol=1e-5, verbose=1):
"""
Check for duplicated entries in the list of ddb files
# Args:
# std_tol: Tolerance on standard deviation
# verbose: Verbosity level.
Args:
std_tol: Tolerance on standard deviation
verbose: Verbosity level.
# Return: (retcode, results) where results maps qpt --> DataFrame with perts as index.
# """
# from pprint import pprint
Return: (retcode, results) where results maps qpt --> DataFrame with perts as index.
"""
from pprint import pprint
# # Build q --> group of dataframes.
# from collections import defaultdict
# q2dfgroup = defaultdict(list)
# for ddb in self.abifiles:
# for qpt, df in ddb.computed_dynmat.items():
# q2dfgroup[qpt].append(df)
# Build q --> group of dataframes.
from collections import defaultdict
q2dfgroup = defaultdict(list)
for ddb in self.abifiles:
for qpt, df in ddb.computed_dynmat.items():
q2dfgroup[qpt].append(df)
# retcode, results = 0, {}
# for qpt, dfgroup in q2dfgroup.items():
# all_indset = [set(df.index) for df in dfgroup]
# # Build union of all dynmat indices with this q
# allps = set(all_indset[0]).union(*all_indset)
# #allps = set(all_indset[0]).intersection(*all_indset)
retcode, results = 0, {}
for qpt, dfgroup in q2dfgroup.items():
all_indset = [set(df.index) for df in dfgroup]
# Build union of all dynmat indices with this q
allps = set(all_indset[0]).union(*all_indset)
#allps = set(all_indset[0]).intersection(*all_indset)
# index, d_list = [], []
# for p in allps:
# # Find dataframes with this p
# found = [p in index for index in all_indset]
# count = found.count(True)
# if count == 1: continue
# if verbose:
# print("Found %s duplicated entries for p: %s" % (count, str(p)))
index, d_list = [], []
for p in allps:
# Find dataframes with this p
found = [p in index for index in all_indset]
count = found.count(True)
if count == 1: continue
if verbose:
print("Found %s duplicated entries for p: %s" % (count, str(p)))
# # Compute stats for this p (complex numbers)
# cvalues = []
# for f, df in zip(found, dfgroup):
# if not f: continue
# c = df["cvalue"].loc[[p]]
# cvalues.append(c)
# Compute stats for this p (complex numbers)
cvalues = []
for f, df in zip(found, dfgroup):
if not f: continue
c = df["cvalue"].loc[[p]]
cvalues.append(c)
# cvalues = np.array(cvalues)
# norms = np.abs(cvalues)
# d = dict(mean=cvalues.mean(), std=cvalues.std(),
# min_norm=norms.min(), max_norm=norms.max(), count=count)
cvalues = np.array(cvalues)
norms = np.abs(cvalues)
d = dict(mean=cvalues.mean(), std=cvalues.std(),
min_norm=norms.min(), max_norm=norms.max(), count=count)
# # Print warning if large deviation
# #if d["max_norm"] - d["min_norm"] > 1e-5:
# if d["std"] > std_tol:
# retcode += 1
# cprint("Found std > %s" % std_tol, "red")
# pprint(cvalues)
# if verbose:
# pprint(d)
# print(2 * "")
# Print warning if large deviation
#if d["max_norm"] - d["min_norm"] > 1e-5:
if d["std"] > std_tol:
retcode += 1
cprint("Found std > %s" % std_tol, "red")
pprint(cvalues)
if verbose:
pprint(d)
print(2 * "")
# d_list.append(d)
# index.append(p)
d_list.append(d)
index.append(p)
# results[qpt] = pd.DataFrame(d_list, index=index)
results[qpt] = pd.DataFrame(d_list, index=index)
return retcode, results
# return retcode, results
def get_dataframe_at_qpoint(self, qpoint=None, units="eV", asr=2, chneut=1, dipdip=1,
with_geo=True, with_spglib=True, abspath=False, funcs=None):
@ -1947,7 +1970,7 @@ class DdbRobot(Robot):
return dict2namedtuple(phbands_plotter=phbands_plotter, phdos_plotter=phdos_plotter)
def anacompare_elastic(self, ddb_header_keys=None, with_structure=True,
with_spglib=True, manager=None, **kwargs):
with_spglib=True, manager=None, verbose=0, **kwargs):
"""
Compute elastic and piezoelectric properties for all DDBs in the robot and build DataFrame.
@ -1957,7 +1980,8 @@ class DdbRobot(Robot):
with_structure: True to add structure parameters to the DataFrame.
with_spglib: True to compute sgplib space group and add it to the DataFrame.
manager: |TaskManager| object. If None, the object is initialized from the configuration file
kwargs: Keyword arguments passed to `ddb.anaget_elastic`.
verbose: verbosity level. Set it to a value > 0 to get more information
kwargs: Keyword arguments passed to `ddb.anaget_elastic`.
Return: DataFrame and list of ElastData objects.
"""
@ -1965,7 +1989,7 @@ class DdbRobot(Robot):
df_list, elastdata_list = [], []
for label, ddb in self.items():
# Invoke anaddb to compute elastic data.
edata = ddb.anaget_elastic(**kwargs)
edata = ddb.anaget_elastic(verbose=verbose, **kwargs)
elastdata_list.append(edata)
# Build daframe with properties derived from the elastic tensor.
@ -1982,20 +2006,83 @@ class DdbRobot(Robot):
# Add path to the DDB file.
df["ddb_path"] = ddb.filepath
df_list.append(df)
# Concatenate dataframes.
return pd.concat(df_list), elastdata_list
def anacompare_becs(self, ddb_header_keys=None, chneut=1, verbose=0):
"""
Compute Born effective charges for all DDBs in the robot and build DataFrame.
with Voigt indices as colums + metadata. Useful for convergence studies.
Args:
ddb_header_keys: List of keywords in the header of the DDB file
whose value will be added to the Dataframe.
verbose: verbosity level. Set it to a value > 0 to get more information
Return: DataFrame and list of Becs objects
"""
ddb_header_keys = [] if ddb_header_keys is None else list_strings(ddb_header_keys)
df_list, becs_list = [], []
for label, ddb in self.items():
# Invoke anaddb to compute Becs
_, becs = ddb.anaget_emacro_and_becs(chneut=chneut, verbose=verbose)
becs_list.append(becs)
df = becs.get_voigt_dataframe()
# Add metadata to the dataframe.
df["chneut"] = chneut
for k in list_strings(ddb_header_keys):
df[k] = ddb.header[k]
# Add path to the DDB file.
df["ddb_path"] = ddb.filepath
df_list.append(df)
# Concatenate dataframes.
return pd.concat(df_list), becs_list
#def anacompare_emacro(self, ddb_header_keys=None, verbose=0):
# """
# Compute Born effective charges for all DDBs in the robot and build DataFrame.
# with Voigt indices as colums + metadata. Useful for convergence studies.
# Args:
# ddb_header_keys: List of keywords in the header of the DDB file
# whose value will be added to the Dataframe.
# verbose: verbosity level. Set it to a value > 0 to get more information
# Return: DataFrame and list of DielectricTensor objects
# """
# ddb_header_keys = [] if ddb_header_keys is None else list_strings(ddb_header_keys)
# df_list, emacro_list = [], []
# for label, ddb in self.items():
# # Invoke anaddb to compute Becs
# emacro, _ = ddb.anaget_emacro_and_becs(chneut=chneut, verbose=verbose)
# emacro_list.append(emacro)
# df = emacro.get_voigt_dataframe()
# # Add metadata to the dataframe.
# for k in list_strings(ddb_header_keys):
# df[k] = ddb.header[k]
# # Add path to the DDB file.
# df["ddb_path"] = ddb.filepath
# df_list.append(df)
# # Concatenate dataframes.
# return pd.concat(df_list), emacro_list
def yield_figs(self, **kwargs): # pragma: no cover
"""
This function *generates* a predefined list of matplotlib figures with minimal input from the user.
"""
print("Invoking anaddb through anaget_phonon_plotters...")
r = self.anaget_phonon_plotters()
for fig in r.phbands_plotter.yield_figs(): yield fig
for fig in r.phdos_plotter.yield_figs(): yield fig
if all(ddb.has_at_least_one_atomic_perturbation() for ddb in self.abifiles):
print("Invoking anaddb through anaget_phonon_plotters...")
r = self.anaget_phonon_plotters()
for fig in r.phbands_plotter.yield_figs(): yield fig
for fig in r.phdos_plotter.yield_figs(): yield fig
def write_notebook(self, nbpath=None):
"""

8
abipy/dfpt/elastic.py
View File

@ -22,9 +22,9 @@ class MyElasticTensor(ElasticTensor):
def _repr_html_(self):
"""Integration with jupyter notebooks."""
return self.to_voigt_dataframe()._repr_html_()
return self.get_voigt_dataframe()._repr_html_()
def to_voigt_dataframe(self, tol=1e-5):
def get_voigt_dataframe(self, tol=1e-5):
tensor = self.zeroed(tol) if tol else self
columns = ["xx", "yy", "zz", "yz", "xz", "xy"]
#columns = ["1", "2", "3", "4", "5", "6"]
@ -39,9 +39,9 @@ class MyPiezoTensor(PiezoTensor):
def _repr_html_(self):
"""Integration with jupyter notebooks."""
return self.to_voigt_dataframe()._repr_html_()
return self.get_voigt_dataframe()._repr_html_()
def to_voigt_dataframe(self, tol=1e-5):
def get_voigt_dataframe(self, tol=1e-5):
tensor = self.zeroed(tol) if tol else self
index = ["Px", "Py", "Pz"]
columns = ["xx", "yy", "zz", "yz", "xz", "xy"]

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

@ -256,6 +256,10 @@ class DdbTest(AbipyTest):
self.assert_almost_equal(emacro, ref_emacro)
repr(becs); str(becs)
assert becs.to_string(verbose=2)
for arr, z in zip(becs.values, becs.zstars):
self.assert_equal(arr, z)
df = becs.get_voigt_dataframe(select_symbols="O")
assert len(df) == 2
# get the dielectric tensor generator from anaddb
dtg = ddb.anaget_dielectric_tensor_generator(verbose=2)
@ -475,14 +479,29 @@ class DdbRobotTest(AbipyTest):
with abilab.DdbRobot.from_files(filepaths) as robot:
robot.add_file("samefile", filepaths[0])
assert len(robot) == 2
# Test anacompare_elastic
ddb_header_keys=["nkpt", "tsmear"]
df, edata_list = robot.anacompare_elastic(ddb_header_keys=ddb_header_keys,
with_structure=True, with_spglib=False, relaxed_ion="automatic", piezo="automatic", verbose=1)
assert "tensor_name" in df.keys()
assert "ddb_path" in df
for k in ddb_header_keys:
assert k in df
assert len(edata_list) == 2
# Test anacompare_becs
df, becs_list = robot.anacompare_becs(ddb_header_keys=ddb_header_keys, chneut=1, verbose=2)
for k in ddb_header_keys:
assert k in df
#assert "ddb_path" in df and (0, 0) in df
assert len(becs_list) == 2
# Test anacompare_emacro
#df, emacro_list = robot.anacompare_emacro(ddb_header_keys=ddb_header_keys, verbose=2)
#assert "ddb_path" in df and (0, 0) in df
#assert len(emacro_list) == 2
class PhononComputationTest(AbipyTest):