Replace old tensor objects with pmg tensors

This commit is contained in:
Matteo Giantomassi 2018-08-08 01:17:26 +02:00
parent 3d8ac10c87
commit a4f35e1302
6 changed files with 211 additions and 265 deletions

View File

@ -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)

View File

@ -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()

View File

@ -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

View File

@ -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)

View File

@ -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.

View File

@ -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
---------------------