Source code for abipy.dfpt.phtk

# coding: utf-8
"""
Phonon Toolkit: This module gathers low-level tools to operate on phonons.
"""
import warnings
import sys
import numpy as np
import abipy.core.abinit_units as abu

from monty.functools import lazy_property
from pymatgen.core.periodic_table import Element
from abipy.core.mixins import Has_Structure
from abipy.iotools import ETSF_Reader


# TODO: amu should become mandatory.
[docs]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 abinit and pymatgen. .. note:: These eigenvectors are orthonormalized and should be very close to the ones computed by Abinit in a.u. Note, however, that the output vectors are given in atomic units so dividing then by the sqrt(Mass) won't give the dipl_cart used in PhononBands that are in Angstrom. Args: phdispl: a numpy array containing the displacements in cartesian coordinates. The last index should have size 3*(num atoms), but the rest of the shape is arbitrary. If qpts is not None the first dimension 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. 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 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_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_symbol[a.specie.symbol]*abu.amu_emass) / abu.Bohr_Ang return eigvec
[docs]def match_eigenvectors(v1, v2): """ Given two list of vectors, returns the pair matching based on the complex scalar product. Returns the indices of the second list that match the vectors of the first list in ascending order. """ prod = np.absolute(np.dot(v1, v2.transpose().conjugate())) indices = np.zeros(len(v1), dtype=np.int) missing_v1 = [True] * len(v1) missing_v2 = [True] * len(v1) for m in reversed(np.argsort(prod, axis=None)): i, j = np.unravel_index(m, prod.shape) if missing_v1[i] and missing_v2[j]: indices[i] = j missing_v1[i] = missing_v2[j] = False if not any(missing_v1): if any(missing_v2): raise RuntimeError('Something went wrong in matching vectors: {} {}'.format(v1, v2)) break return indices
[docs]class NonAnalyticalPh(Has_Structure): """ Phonon data at gamma including non analytical contributions Read from anaddb.nc """ def __init__(self, structure, directions, phfreqs, phdispl_cart, amu=None): """ Args: structure: |Structure| object. directions: Cartesian directions along which the non analytical frequencies have been calculated phfreqs: Phonon frequencies with non analytical contribution in eV along directions phdispl_cart: Displacement in Angstrom in Cartesian coordinates with non analytical contribution along directions amu: dictionary that associates the atomic species present in the structure to the values of the atomic mass units used for the calculation """ self._structure = structure self.directions = directions self.phfreqs = phfreqs self.phdispl_cart = phdispl_cart self.amu = amu self.amu_symbol = None if amu is not None: self.amu_symbol = {} for z, m in amu.items(): el = Element.from_Z(int(z)) self.amu_symbol[el.symbol] = m
[docs] @classmethod def from_file(cls, filepath): """ Reads the non analytical directions, frequencies and displacements from the anaddb.nc file specified. Non existence of displacements is accepted for compatibility with abinit 8.0.6 Raises an error if the other values are not present in anaddb.nc. """ with ETSF_Reader(filepath) as r: directions = r.read_value("non_analytical_directions") phfreq = r.read_value("non_analytical_phonon_modes") # need a default as the first abinit version including IFCs in the netcdf doesn't have this attribute phdispl_cart = r.read_value("non_analytical_phdispl_cart", cmode="c", default=None) structure = r.read_structure() amu_list = r.read_value("atomic_mass_units", default=None) if amu_list is not None: # ntypat arrays atomic_numbers = r.read_value("atomic_numbers") amu = {at: a for at, a in zip(atomic_numbers, amu_list)} else: amu = None return cls(structure=structure, directions=directions, phfreqs=phfreq, phdispl_cart=phdispl_cart, amu=amu)
[docs] @lazy_property def dyn_mat_eigenvect(self): """ [ndirection, 3*natom, 3*natom] array with the orthonormal eigenvectors of the dynamical matrix. in Cartesian coordinates. """ return get_dyn_mat_eigenvec(self.phdispl_cart, self.structure, amu=self.amu)
@property def structure(self): """|Structure| object.""" return self._structure
[docs] def index_direction(self, direction, cartesian=False): """ Returns: the index of direction. Raises: `ValueError` if not found. Args: direction: a 3 element list indicating the direction. Can be a generic vector cartesian: if True the direction are already in cartesian coordinates, if False it will be converted to match the internal description of the directions. """ if not cartesian: direction = self.structure.lattice.reciprocal_lattice_crystallographic.get_cartesian_coords(direction) else: direction = np.array(direction) direction = direction / np.linalg.norm(direction) for i, d in enumerate(self.directions): d = d / np.linalg.norm(d) if np.allclose(d, direction): return i raise ValueError("Cannot find direction: `%s` with cartesian: `%s` in non_analytical cartesian directions:\n%s" % (str(direction), cartesian, str(self.directions)))
[docs] def has_direction(self, direction, cartesian=False): """ Checks if the input direction is among those available. Args: direction: a 3 element list indicating the direction. Can be a generic vector cartesian: if True the direction are already in cartesian coordinates, if False it will be converted to match the internal description of the directions. """ try: self.index_direction(direction, cartesian=cartesian) return True except ValueError: return False
[docs]def open_file_phononwebsite(filename, port=8000, website="http://henriquemiranda.github.io/phononwebsite", host="localhost", browser=None): # pragma: no cover """ Take a file, detect the type and open it on the phonon website Based on a similar function in <https://github.com/henriquemiranda/phononwebsite/phononweb.py> Args: filename: file with phonon data in phononwebsite format. port: Initial port. website: Website URL host: localhost name. browser: Open webpage in ``browser``. Use default if $BROWSER if None. """ if filename.endswith(".json"): filetype = "json" elif filename.endswith(".yaml"): filetype = "yaml" else: filetype = "rest" from http.server import HTTPServer, SimpleHTTPRequestHandler # Add CORS header to the website class CORSRequestHandler(SimpleHTTPRequestHandler): def end_headers(self): #self.send_header('Access-Control-Allow-Origin', website) self.send_header('Access-Control-Allow-Origin', "http://henriquemiranda.github.io") SimpleHTTPRequestHandler.end_headers(self) def log_message(self, format, *args): return # Initialize http server thread print('Starting HTTP server at port %d ...' % port, end=" ") trial, max_ntrial = 0, 50 while trial < max_ntrial: try: server = HTTPServer(('', port), CORSRequestHandler) #print("got port:", port) break except OSError: trial += 1 port += 1 print(port, end=", ") else: raise RuntimeError("Cannot find available port after %s attempts" % max_ntrial) # Create threads python server.url = 'http://{}:{}'.format(host, server.server_port) from threading import Thread t = Thread(target=server.serve_forever) t.daemon = True t.start() # Open website with the file try: from urllib.parse import quote except ImportError: from urllib import quote url_filename = 'http://{}:{}/{}'.format(host, server.server_port, quote(filename)) url = '%s/phonon.html?%s=%s' % (website, filetype, url_filename) print("\nOpening URL:", url) print("Using default browser, if the webpage is not displayed correctly", "\ntry to change browser either via command line options or directly in the shell with e.g:\n\n" " export BROWSER=firefox\n") print('Press Ctrl+C to terminate the HTTP server') import webbrowser webbrowser.get(browser).open_new_tab(url) def signal_handler(signal, frame): """Quit application when SIGINT is received""" sys.exit(0) import signal signal.signal(signal.SIGINT, signal_handler) signal.pause()