diff --git a/abipy/core/tensor.py b/abipy/core/tensor.py deleted file mode 100644 index 95573095..00000000 --- a/abipy/core/tensor.py +++ /dev/null @@ -1,173 +0,0 @@ -# coding: utf-8 -""" -This module contains classes representing tensors and helper functions to change lattice. -""" -from __future__ import print_function, division, unicode_literals, absolute_import - -import itertools -import numpy as np - -from monty.dev import deprecated -from pymatgen.symmetry.analyzer import SpacegroupAnalyzer -from abipy.core import Structure - - -__all__ = [ - "Tensor", - "SymmetricTensor", -] - - -def from_cart_to_red(cartesian_tensor,lattice): - mat = lattice.inv_matrix - red_tensor = np.dot(np.dot(np.transpose(mat), cartesian_tensor), mat) - return red_tensor - - -class Tensor(object): - """Representation of a 3x3 tensor""" - - # TODO Remove - #@deprecated(message="abipy.core.Tensor is deprecated and will be replaced by pymatgen tensor in v0.4") - def __init__(self, red_tensor, lattice, space="r"): - """ - Args: - red_tensor: array-like object with the 9 cartesian components of the tensor - lattice: Lattice object defining the reference system - space: - "r" if the lattice is a real space lattice - "g" if the lattice is a reciprocal space lattice - """ - self._reduced_tensor = red_tensor - self._lattice = lattice - self.space = space - - if space == "g": - self._is_real_space = False - elif space == "r": - self._is_real_space = True - else: - raise ValueError("space should be either 'g' or 'r'") - - def __eq__(self, other): - if other is None: return False - return (np.allclose(self.reduced_tensor, other.reduced_tensor) and - self.lattice == other.lattice and - self.space == other.space) - - def __ne__(self, other): - return not (self == other) - - def __repr__(self): - return self.to_string() - - def __str__(self): - return repr(self) - - def to_string(self, verbose=0, with_reduced=False): - lines = [] - app = lines.append - - app("Tensor in %s space." % self.space) - app("") - app("Cartesian coordinates:") - app(str(self.cartesian_tensor)) - - if with_reduced: - app("") - app(str(self.lattice)) - app("Reduced coordinates:") - app(str(self.reduced_tensor)) - - return "\n".join(lines) - - @property - def lattice(self): - return self._lattice - - @property - def reduced_tensor(self): - return self._reduced_tensor - - @property - def is_real_space(self): - return self._is_real_space - - @property - def cartesian_tensor(self): - mat = self._lattice.matrix - return np.dot(np.dot(np.transpose(mat), self._reduced_tensor), mat) - - @classmethod - def from_cartesian_tensor(cls, cartesian_tensor, lattice, space="r"): - red_tensor = from_cart_to_red(cartesian_tensor, lattice) - return cls(red_tensor, lattice,space) - - def symmetrize(self, structure): - tensor = self._reduced_tensor - - if self._is_real_space: - real_lattice = self._lattice - else: - real_lattice = self._lattice.reciprocal_lattice - - # I guess this is the reason why tensor.symmetrize (omega) is so slow! - real_finder = SpacegroupAnalyzer(structure) - - real_symmops = real_finder.get_point_group_operations(cartesian=True) - - cartesian_tensor = self.cartesian_tensor - - sym_tensor = np.zeros((3,3)) - - my_tensor = cartesian_tensor - - for real_sym in real_symmops: - mat = real_sym.rotation_matrix - prod_sym = np.dot(np.transpose(mat),np.dot(cartesian_tensor,mat)) - sym_tensor = sym_tensor + prod_sym - - sym_tensor = sym_tensor/len(real_symmops) - - self._reduced_tensor = from_cart_to_red(sym_tensor,self._lattice) - - -class SymmetricTensor(Tensor): - """Representation of a 3x3 symmetric tensor""" - - @classmethod - def from_directions(cls, qpoints, values, lattice, space): - """ - Build a `SymmetricTensor` from the values computed along 6 directions. - - Args: - qpoints: fractional coordinates of 6 independent q-directions - values: values of (q^T E q)/(q^T q) along the 6 qpoints - lattice: `Lattice` object defining the reference system - space: "r" if the lattice is a real space lattice - "g" if the lattice is a reciprocal space lattice - """ - assert len(qpoints) == 6 and len(values) == len(qpoints) - - mat = lattice.matrix - metric = np.dot(np.transpose(mat),mat) - - coeffs_red = np.zeros((6,6)) - - for (iqpt,qpt) in enumerate(qpoints): - metqpt = np.dot(metric,qpt) - - coeffs_red[iqpt,:] = [metqpt[0]**2,metqpt[1]**2,metqpt[2]**2, - 2*metqpt[0]*metqpt[1],2*metqpt[0]*metqpt[2],2*metqpt[1]*metqpt[2]] - - normqpt_red = np.dot(np.transpose(qpt),np.dot(metric,qpt)) - - coeffs_red[iqpt,:] = coeffs_red[iqpt,:] / normqpt_red - - red_symm = np.linalg.solve(coeffs_red,values) - - red_tensor = [[red_symm[0],red_symm[3],red_symm[4]], - [red_symm[3],red_symm[1],red_symm[5]], - [red_symm[4],red_symm[5],red_symm[2]]] - - return cls(red_tensor, lattice, space) diff --git a/abipy/core/tests/test_tensor.py b/abipy/core/tests/test_tensor.py deleted file mode 100644 index 234044ea..00000000 --- a/abipy/core/tests/test_tensor.py +++ /dev/null @@ -1,52 +0,0 @@ -"""Tests for tensor""" -from __future__ import print_function, division - -import numpy as np - -from pymatgen.core.structure import Structure -from pymatgen.core.lattice import Lattice -from pymatgen.symmetry.analyzer import SpacegroupAnalyzer -from abipy.core.tensor import * -from abipy.core.testing import AbipyTest - - -class TestTensor(AbipyTest): - """Unit tests for Tensor.""" - - def test_tensor(self): - """Initialize Tensor""" - - lattice = Lattice.hexagonal(4,6) - #rprimd = np.array([[0,0.5,0.5],[0.5,0,0.5],[0.5,0.5,0]]) - #rprimd = rprimd*10 - #lattice = Lattice(rprimd) - structure = Structure(lattice, ["Ga", "As"], - [[0, 0, 0], [0.5, 0.5, 0.5]]) - - #finder = SymmetryFinder(structure) - finder = SpacegroupAnalyzer(structure) - - spacegroup = finder.get_space_group_operations() - pointgroup = finder.get_point_group_symbol() - - cartesian_tensor = [[2,3,1.2],[3,4,1.0],[1.2,1.0,6]] - - tensor = Tensor.from_cartesian_tensor(cartesian_tensor,lattice.reciprocal_lattice,space="g") - red_tensor = tensor.reduced_tensor - tensor2 = Tensor(red_tensor,lattice.reciprocal_lattice,space="g") - assert(((np.abs(tensor2.cartesian_tensor)-np.abs(cartesian_tensor)) < 1E-8).all()) - - self.assertTrue(tensor==tensor2) - print(tensor) - - #print("non-symmetrized cartesian_tensor = ",tensor2.cartesian_tensor) - tensor2.symmetrize(structure) - - #print("symmetrized_cartesian_tensor = ",tensor2.cartesian_tensor) - - self.serialize_with_pickle(tensor) - - -if __name__ == "__main__": - import unittest - unittest.main() diff --git a/abipy/dfpt/anaddbnc.py b/abipy/dfpt/anaddbnc.py index e27cd7b5..f749b085 100644 --- a/abipy/dfpt/anaddbnc.py +++ b/abipy/dfpt/anaddbnc.py @@ -11,14 +11,13 @@ from collections import OrderedDict from monty.functools import lazy_property from monty.string import marquee from monty.termcolor import cprint -from abipy.core.tensor import Tensor from abipy.core.mixins import AbinitNcFile, Has_Structure, NotebookWriter from abipy.abio.robots import Robot from abipy.iotools import ETSF_Reader from abipy.tools.plotting import add_fig_kwargs, get_axarray_fig_plt, rotate_ticklabels, set_visible from abipy.dfpt.phonons import InteratomicForceConstants from abipy.dfpt.ddb import Becs -from abipy.dfpt.tensors import NLOpticalSusceptibilityTensor +from abipy.dfpt.tensors import Tensor, DielectricTensor, NLOpticalSusceptibilityTensor from abipy.dfpt.elastic import ElasticData @@ -91,6 +90,11 @@ class AnaddbNcFile(AbinitNcFile, Has_Structure, NotebookWriter): app("") app(self.structure.to_string(verbose=verbose, title="Structure")) + import json + app(marquee("Parameters", mark="=")) + app(json.dumps(self.params, indent=2, sort_keys=True)) + app("") + if self.has_elastic_data: app("") df = self.elastic_data.get_average_elastic_dataframe(tensor="elastic_relaxed") @@ -104,38 +108,58 @@ class AnaddbNcFile(AbinitNcFile, Has_Structure, NotebookWriter): app(df.T.to_string(index=True)) if verbose: - df = self.elastic_data.get_voigt_dataframe(["elastic_relaxed", "elastic_clamped"]) + df = self.elastic_data.get_voigt_dataframe() app(df.T.to_string()) + if self.emacro is not None: + app("Macroscopic dielectric tensor (Cartesian coords)") + app(str(self.emacro)) + app("") + + if self.emacro_rlx is not None: + app("Relaxed ion Macroscopic dielectric tensor (Cartesian coords)") + app(str(self.emacro_rlx)) + app("") + + #if self.becs is not None: + + if self.dchide is not None: + app("Non-linear optical susceptibility tensor.") + app(str(self.dchide)) + app("") + + if self.dchidt is not None: + app("First-order change in the linear dielectric susceptibility.") + app(str(self.dchidt)) + app("") + #if self.has_piezoelectric_data: # df = self.elastic_data.get_piezoelectric_dataframe() - #app(str(self.params)) - return "\n".join(lines) @lazy_property def emacro(self): """ - Macroscopic dielectric tensor. None if the file does not contain this information. + Macroscopic dielectric tensor. + None if the file does not contain this information. """ try: - return Tensor.from_cartesian_tensor(self.reader.read_value("emacro_cart"), - self.structure.lattice, space="r") + return DielectricTensor(self.reader.read_value("emacro_cart")) except Exception as exc: - print(exc, "Returning None", sep="\n") + #print(exc, "Returning None", sep="\n") return None @lazy_property def emacro_rlx(self): """ - Relaxed ion Macroscopic dielectric tensor. None if the file does not contain this information. + Relaxed ion Macroscopic dielectric tensor. + None if the file does not contain this information. """ try: - return Tensor.from_cartesian_tensor(self.reader.read_value("emacro_cart_rlx"), - self.structure.lattice, space="r") + return DielectricTensor(self.reader.read_value("emacro_cart_rlx")) except Exception as exc: - print(exc, "Requires dieflag > 0", "Returning None", sep="\n") + #print(exc, "Requires dieflag > 0", "Returning None", sep="\n") return None @lazy_property @@ -147,7 +171,7 @@ class AnaddbNcFile(AbinitNcFile, Has_Structure, NotebookWriter): try: return Becs(self.reader.read_value("becs_cart"), self.structure, chneut=chneut, order="f") except Exception as exc: - print(exc, "Returning None", sep="\n") + #print(exc, "Returning None", sep="\n") return None @lazy_property @@ -160,8 +184,8 @@ class AnaddbNcFile(AbinitNcFile, Has_Structure, NotebookWriter): try: return InteratomicForceConstants.from_file(self.filepath) except Exception as exc: - print(exc) - cprint("Interatomic force constants have not been calculated. Returning None", "red") + #print(exc) + #cprint("Interatomic force constants have not been calculated. Returning None", "red") return None @lazy_property @@ -173,7 +197,7 @@ class AnaddbNcFile(AbinitNcFile, Has_Structure, NotebookWriter): try: return NLOpticalSusceptibilityTensor(self.reader.read_value("dchide")) except Exception as exc: - print(exc, "Requires nlflag > 0", "Returning None", sep="\n") + #print(exc, "Requires nlflag > 0", "Returning None", sep="\n") return None @lazy_property @@ -186,19 +210,20 @@ class AnaddbNcFile(AbinitNcFile, Has_Structure, NotebookWriter): None if the file does not contain this information. """ try: - a = self.reader.read_value("dchidt").T - dchidt = [] - for i in a: - d = [] - for j in i: - d.append(Tensor.from_cartesian_tensor(j, self.structure.lattice, space="r")) - dchidt.append(d) - - return dchidt + a = self.reader.read_value("dchidt").T.copy() except Exception as exc: - print(exc, "Requires 0 < nlflag < 3", "Returning None", sep="\n") + #print(exc, "Requires 0 < nlflag < 3", "Returning None", sep="\n") return None + dchidt = [] + for i in a: + d = [] + for j in i: + d.append(Tensor(j)) + dchidt.append(d) + + return dchidt + @lazy_property def oscillator_strength(self): """ @@ -209,7 +234,7 @@ class AnaddbNcFile(AbinitNcFile, Has_Structure, NotebookWriter): try: return self.reader.read_value("oscillator_strength", cmode="c") except Exception as exc: - print(exc, "Oscillator strengths require dieflag == 1, 3 or 4", "Returning None", sep="\n") + #print(exc, "Oscillator strengths require dieflag == 1, 3 or 4", "Returning None", sep="\n") return None @lazy_property @@ -313,7 +338,6 @@ class AnaddbNcRobot(Robot): rows.append(d) - #index = index if not abspath else self._to_relpaths(index) return pd.DataFrame(rows, index=index, columns=list(rows[0].keys() if rows else None)) @add_fig_kwargs diff --git a/abipy/electrons/bse.py b/abipy/electrons/bse.py index c37cf675..ee85aa51 100644 --- a/abipy/electrons/bse.py +++ b/abipy/electrons/bse.py @@ -16,7 +16,6 @@ from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, get_axarray_fig from abipy.core.func1d import Function1D from abipy.core.kpoints import Kpoint, KpointList from abipy.core.mixins import AbinitNcFile, Has_Structure, NotebookWriter -from abipy.core.tensor import SymmetricTensor from abipy.iotools import ETSF_Reader from abipy.tools.plotting import set_axlims from abipy.tools import duck @@ -54,7 +53,7 @@ class DielectricTensor(object): # One tensor for each frequency all_tensors = [] for ifrq, freq in enumerate(mdf.wmesh): - tensor = SymmetricTensor.from_directions(mdf.qfrac_coords, all_emacros[:,ifrq], + tensor = _SymmetricTensor.from_directions(mdf.qfrac_coords, all_emacros[:,ifrq], structure.lattice.reciprocal_lattice, space="g") all_tensors.append(tensor) @@ -1023,3 +1022,159 @@ class MdfRobot(Robot, RobotWithEbands): nb.cells.extend(self.get_ebands_code_cells()) return self._write_nb_nbpath(nb, nbpath) + + +def _from_cart_to_red(cartesian_tensor,lattice): + mat = lattice.inv_matrix + red_tensor = np.dot(np.dot(np.transpose(mat), cartesian_tensor), mat) + return red_tensor + + +# TODO Remove +#@deprecated(message="abipy.core.Tensor is deprecated and will be replaced by pymatgen tensor in v0.4") +class _Tensor(object): + """Representation of a 3x3 tensor""" + + def __init__(self, red_tensor, lattice, space="r"): + """ + Args: + red_tensor: array-like object with the 9 cartesian components of the tensor + lattice: Lattice object defining the reference system + space: + "r" if the lattice is a real space lattice + "g" if the lattice is a reciprocal space lattice + """ + self._reduced_tensor = red_tensor + self._lattice = lattice + self.space = space + + if space == "g": + self._is_real_space = False + elif space == "r": + self._is_real_space = True + else: + raise ValueError("space should be either 'g' or 'r'") + + def __eq__(self, other): + if other is None: return False + return (np.allclose(self.reduced_tensor, other.reduced_tensor) and + self.lattice == other.lattice and + self.space == other.space) + + def __ne__(self, other): + return not (self == other) + + def __repr__(self): + return self.to_string() + + def __str__(self): + return repr(self) + + def to_string(self, verbose=0, with_reduced=False): + lines = [] + app = lines.append + + app("Tensor in %s space." % self.space) + app("") + app("Cartesian coordinates:") + app(str(self.cartesian_tensor)) + + if with_reduced: + app("") + app(str(self.lattice)) + app("Reduced coordinates:") + app(str(self.reduced_tensor)) + + return "\n".join(lines) + + @property + def lattice(self): + return self._lattice + + @property + def reduced_tensor(self): + return self._reduced_tensor + + @property + def is_real_space(self): + return self._is_real_space + + @property + def cartesian_tensor(self): + mat = self._lattice.matrix + return np.dot(np.dot(np.transpose(mat), self._reduced_tensor), mat) + + @classmethod + def from_cartesian_tensor(cls, cartesian_tensor, lattice, space="r"): + red_tensor = _from_cart_to_red(cartesian_tensor, lattice) + return cls(red_tensor, lattice,space) + + def symmetrize(self, structure): + tensor = self._reduced_tensor + + if self._is_real_space: + real_lattice = self._lattice + else: + real_lattice = self._lattice.reciprocal_lattice + + # I guess this is the reason why tensor.symmetrize (omega) is so slow! + from pymatgen.symmetry.analyzer import SpacegroupAnalyzer + real_finder = SpacegroupAnalyzer(structure) + + real_symmops = real_finder.get_point_group_operations(cartesian=True) + + cartesian_tensor = self.cartesian_tensor + + sym_tensor = np.zeros((3,3)) + + my_tensor = cartesian_tensor + + for real_sym in real_symmops: + mat = real_sym.rotation_matrix + prod_sym = np.dot(np.transpose(mat),np.dot(cartesian_tensor,mat)) + sym_tensor = sym_tensor + prod_sym + + sym_tensor = sym_tensor/len(real_symmops) + + self._reduced_tensor = _from_cart_to_red(sym_tensor,self._lattice) + + +class _SymmetricTensor(_Tensor): + """Representation of a 3x3 symmetric tensor""" + + @classmethod + def from_directions(cls, qpoints, values, lattice, space): + """ + Build a `_SymmetricTensor` from the values computed along 6 directions. + + Args: + qpoints: fractional coordinates of 6 independent q-directions + values: values of (q^T E q)/(q^T q) along the 6 qpoints + lattice: `Lattice` object defining the reference system + space: "r" if the lattice is a real space lattice + "g" if the lattice is a reciprocal space lattice + """ + assert len(qpoints) == 6 and len(values) == len(qpoints) + + mat = lattice.matrix + metric = np.dot(np.transpose(mat),mat) + + coeffs_red = np.zeros((6,6)) + + for (iqpt,qpt) in enumerate(qpoints): + metqpt = np.dot(metric,qpt) + + coeffs_red[iqpt,:] = [metqpt[0]**2,metqpt[1]**2,metqpt[2]**2, + 2*metqpt[0]*metqpt[1],2*metqpt[0]*metqpt[2],2*metqpt[1]*metqpt[2]] + + normqpt_red = np.dot(np.transpose(qpt),np.dot(metric,qpt)) + + coeffs_red[iqpt,:] = coeffs_red[iqpt,:] / normqpt_red + + red_symm = np.linalg.solve(coeffs_red,values) + + red_tensor = [[red_symm[0],red_symm[3],red_symm[4]], + [red_symm[3],red_symm[1],red_symm[5]], + [red_symm[4],red_symm[5],red_symm[2]]] + + return cls(red_tensor, lattice, space) diff --git a/abipy/integration_tests/TODO.md b/abipy/integration_tests/TODO.md index 1d3ce0a4..fb0834d2 100644 --- a/abipy/integration_tests/TODO.md +++ b/abipy/integration_tests/TODO.md @@ -74,7 +74,7 @@ TODO list: * Scheduler should report info on exceptions (especially if at the end when on_all_ok is invoked) -* Replace core.tensor with pymatgen tensor (postponed to v0.4) +* ALMOST DONE: Replace core.tensor with pymatgen tensor Check DielectricTensor in Anaddb from DDB. Use pmg tensor for stress as well. diff --git a/docs/api/core_api.rst b/docs/api/core_api.rst index b7f9ee7e..d55fefac 100644 --- a/docs/api/core_api.rst +++ b/docs/api/core_api.rst @@ -116,14 +116,6 @@ core Package :undoc-members: :show-inheritance: -:mod:`tensor` Module --------------------- - -.. automodule:: abipy.core.tensor - :members: - :undoc-members: - :show-inheritance: - :mod:`testing` Module ---------------------