# Path: pyDFTutils/perovskite/perovskite_mode.py
# def label_zone_boundary(qname,
#                         phdisp=None,
#                         masses=None,
#                         evec=None,
#                         notation='IR'):
#     IR_translation = {}
# 
#     IR_translation['Gamma'] = {
#         '$\Delta_1$': r'$\Gamma_4^-$',
#         '$\Delta_2$': r'',
#         '$\Delta_5$': r'',
#     }
# 
#     IR_translation['R'] = {
#         r'$\Gamma_2\prime$': '$R_2^-$',
#         r'$\Gamma_{12}\prime$': '$R_3^-$',
#         r'$\Gamma_{25}$': '$R_5^-$',
#         r'$\Gamma_{25}\prime$': '$R_5^+$',
#         r'$\Gamma_{15}$': '$R_4^-$',
#     }
# 
#     IR_translation['X'] = {
#         '$M_1$': '$X_1^+$',
#         '$M_2\prime$': '$X_3^-$',
#         '$M_3$': '$X_2^+$',
#         '$M_5$': '$X_5^+$',
#         '$M_5\prime$': '$X_5^-$',
#     }
# 
#     IR_translation['M'] = {
#         '$M_1$': '$M_1^+$',
#         '$M_2$': '$M_3^+$',
#         '$M_3$': '$M_2^+$',
#         '$M_4$': '$M_4^+$',
#         '$M_2\prime$': '$M_3^-$',
#         '$M_3\prime$': '$M_2^-$',
#         '$M_5$': '$M_5^+$',
#         '$M_5\prime$': '$M_5^-$',
#     }
# 
#     if phdisp is not None or masses is not None:
#         evec = np.array(phdisp) * np.sqrt(np.kron(masses, [1, 1, 1]))
#         evec = np.real(evec) / np.linalg.norm(evec)
# 
#     mode = None
#     for m in IR_dict[qname]:
#         mvec = np.real(m)
#         mvec = mvec / np.linalg.norm(mvec)
#         evec = evec / np.linalg.norm(evec)
#         p = np.abs(np.dot(np.real(evec), mvec))
#         if p > 0.4:  #1.0 / np.sqrt(2):
#             print("-------------")
#             #print("Found! p= %s" % p)
#             #print("eigen vector: ", nmode._make(mvec))
#             if notation.lower() == 'cowley':
#                 mode = IR_dict[qname][m]
#             elif notation.upper() == 'IR':
#                 #print(IR_translation[qname])
#                 mode = IR_translation[qname][IR_dict[qname][m]]
#             elif notation.lower() == 'both':
#                 mode = (IR_dict[qname][m],
#                         IR_translation[qname][IR_dict[qname][m]])
#             else:
#                 raise ValueError('notation should be Cowley|IR|both')
#             #print("mode: ", mode, m)
#     if mode is None:
#         print("==============")
#         print("eigen vector: ", nmode._make(evec))
#     return mode
# 
# def label_Gamma(phdisp=None, masses=None, evec=None):
#     if phdisp is not None and masses is not None:
#         evec = np.array(phdisp) * np.sqrt(np.kron(masses, [1, 1, 1]))
#         evec = np.real(evec) / np.linalg.norm(evec)
# 
#     basis_dict = mass_to_Gamma_basis_3d(masses, eigen_type='eigen_vector')
#     mode = None
# 
#     for m in basis_dict:
#         mvec = np.real(m)
#         mvec = mvec / np.linalg.norm(mvec)
#         evec = evec / np.linalg.norm(evec)
#         p = np.abs(np.dot(np.real(evec), mvec))
#         if p > 0.5:  #1.0 / np.sqrt(2):
#             mode = basis_dict[m]
#     if mode is None:
#         print("==============")
#         print("eigen vector: ", nmode._make(evec))
#     return mode


import os
import numpy as np
import matplotlib.pyplot as plt
from ase.data import chemical_symbols
from abipy.abilab import abiopen
from pyDFTutils.perovskite.perovskite_mode import label_zone_boundary, label_Gamma
from ase.units import Ha
from spglib import spglib
    def get_gamma_modes(self):
        """
        return (Freqs, names, evecs)
        """
        return self.phonon_mode_freqs['Gamma'], self.phonon_mode_names['Gamma'], self.phonon_mode_evecs['Gamma'], 


    def get_gamma_mode(self, mode_name):
        """
        return the frequencies of mode name.
        """
        ibranches = []
        freqs = []
        for imode, mode in enumerate(zip(self.phonon_mode_freqs['Gamma'], self.phonon_mode_names['Gamma'])):
            freq, mname = mode
            if mname == mode_name:
                ibranches.append(imode)
                freqs.append(freq)
        return ibranches, freqs

    def label_Gamma_all(self, qpoints, evals, evecs, label=True):
        Gamma_mode_freqs = []
        Gamma_mode_names = []
        Gamma_mode_evecs = []
        for i, qpt in enumerate(qpoints):
            if np.isclose(qpt, [0, 0, 0], rtol=1e-5, atol=1e-3).all():
                evecq = evecs[i]
                for j, evec in enumerate(evecq.T):
                    freq = evals[i][j]
                    if label:

# Path: pyDFTutils/perovskite/perovskite_mode.py
# def label_zone_boundary(qname,
#                         phdisp=None,
#                         masses=None,
#                         evec=None,
#                         notation='IR'):
#     IR_translation = {}
# 
#     IR_translation['Gamma'] = {
#         '$\Delta_1$': r'$\Gamma_4^-$',
#         '$\Delta_2$': r'',
#         '$\Delta_5$': r'',
#     }
# 
#     IR_translation['R'] = {
#         r'$\Gamma_2\prime$': '$R_2^-$',
#         r'$\Gamma_{12}\prime$': '$R_3^-$',
#         r'$\Gamma_{25}$': '$R_5^-$',
#         r'$\Gamma_{25}\prime$': '$R_5^+$',
#         r'$\Gamma_{15}$': '$R_4^-$',
#     }
# 
#     IR_translation['X'] = {
#         '$M_1$': '$X_1^+$',
#         '$M_2\prime$': '$X_3^-$',
#         '$M_3$': '$X_2^+$',
#         '$M_5$': '$X_5^+$',
#         '$M_5\prime$': '$X_5^-$',
#     }
# 
#     IR_translation['M'] = {
#         '$M_1$': '$M_1^+$',
#         '$M_2$': '$M_3^+$',
#         '$M_3$': '$M_2^+$',
#         '$M_4$': '$M_4^+$',
#         '$M_2\prime$': '$M_3^-$',
#         '$M_3\prime$': '$M_2^-$',
#         '$M_5$': '$M_5^+$',
#         '$M_5\prime$': '$M_5^-$',
#     }
# 
#     if phdisp is not None or masses is not None:
#         evec = np.array(phdisp) * np.sqrt(np.kron(masses, [1, 1, 1]))
#         evec = np.real(evec) / np.linalg.norm(evec)
# 
#     mode = None
#     for m in IR_dict[qname]:
#         mvec = np.real(m)
#         mvec = mvec / np.linalg.norm(mvec)
#         evec = evec / np.linalg.norm(evec)
#         p = np.abs(np.dot(np.real(evec), mvec))
#         if p > 0.4:  #1.0 / np.sqrt(2):
#             print("-------------")
#             #print("Found! p= %s" % p)
#             #print("eigen vector: ", nmode._make(mvec))
#             if notation.lower() == 'cowley':
#                 mode = IR_dict[qname][m]
#             elif notation.upper() == 'IR':
#                 #print(IR_translation[qname])
#                 mode = IR_translation[qname][IR_dict[qname][m]]
#             elif notation.lower() == 'both':
#                 mode = (IR_dict[qname][m],
#                         IR_translation[qname][IR_dict[qname][m]])
#             else:
#                 raise ValueError('notation should be Cowley|IR|both')
#             #print("mode: ", mode, m)
#     if mode is None:
#         print("==============")
#         print("eigen vector: ", nmode._make(evec))
#     return mode
# 
# def label_Gamma(phdisp=None, masses=None, evec=None):
#     if phdisp is not None and masses is not None:
#         evec = np.array(phdisp) * np.sqrt(np.kron(masses, [1, 1, 1]))
#         evec = np.real(evec) / np.linalg.norm(evec)
# 
#     basis_dict = mass_to_Gamma_basis_3d(masses, eigen_type='eigen_vector')
#     mode = None
# 
#     for m in basis_dict:
#         mvec = np.real(m)
#         mvec = mvec / np.linalg.norm(mvec)
#         evec = evec / np.linalg.norm(evec)
#         p = np.abs(np.dot(np.real(evec), mvec))
#         if p > 0.5:  #1.0 / np.sqrt(2):
#             mode = basis_dict[m]
#     if mode is None:
#         print("==============")
#         print("eigen vector: ", nmode._make(evec))
#     return mode


# Path: pyDFTutils/phonon/parser.py
import os
import numpy as np
import matplotlib.pyplot as plt
from ase.data import chemical_symbols
from abipy.abilab import abiopen
from pyDFTutils.perovskite.perovskite_mode import label_zone_boundary, label_Gamma
from ase.units import Ha
from spglib import spglib

    def get_gamma_modes(self):
        """
        return (Freqs, names, evecs)
        """
        return self.phonon_mode_freqs['Gamma'], self.phonon_mode_names['Gamma'], self.phonon_mode_evecs['Gamma'], 


    def get_gamma_mode(self, mode_name):
        """
        return the frequencies of mode name.
        """
        ibranches = []
        freqs = []
        for imode, mode in enumerate(zip(self.phonon_mode_freqs['Gamma'], self.phonon_mode_names['Gamma'])):
            freq, mname = mode
            if mname == mode_name:
                ibranches.append(imode)
                freqs.append(freq)
        return ibranches, freqs

    def label_Gamma_all(self, qpoints, evals, evecs, label=True):
        Gamma_mode_freqs = []
        Gamma_mode_names = []
        Gamma_mode_evecs = []
        for i, qpt in enumerate(qpoints):
            if np.isclose(qpt, [0, 0, 0], rtol=1e-5, atol=1e-3).all():
                evecq = evecs[i]
                for j, evec in enumerate(evecq.T):
                    freq = evals[i][j]
                    if label:

                        mode = label_Gamma(