"""
Useful functions used in embedding_ifc
"""
from __future__ import annotations
import os
import shutil
import numpy as np
from abipy.core.structure import Structure
[docs]
def stru_0_1_to_minus_05_05(structure: Structure):
"""
Translate the structure so that the frac_coords go from [0,1] to [-0.5,0,5]
"""
new_stru = structure.copy()
for site in new_stru:
if site.frac_coords[0] >= 0.5:
site.frac_coords[0] = site.frac_coords[0] - 1
if site.frac_coords[1] >= 0.5:
site.frac_coords[1] = site.frac_coords[1] - 1
if site.frac_coords[2] >= 0.5:
site.frac_coords[2] = site.frac_coords[2] - 1
return new_stru
[docs]
def frac_coords_05_to_minus05(structure: Structure) -> Structure:
"""
Move atoms that are close to frac_coord = 0.5 to -0.5
Needed for atoms that are at -0.5 in the perfect supercell but move to for instance 0.498 after relaxation
"""
new_stru = structure.copy()
for i, atom in enumerate(new_stru):
if atom.frac_coords[0] > 0.495:
atom.frac_coords[0] = -0.5
if atom.frac_coords[1] > 0.495:
atom.frac_coords[1] = -0.5
if atom.frac_coords[2] > 0.495:
atom.frac_coords[2] = -0.5
return new_stru
[docs]
def center_wrt_defect(structure: Structure, defect_coord) -> Structure:
"""
Center the structure around defect_coord. Defect is now at [0,0,0]
"""
new_stru = structure.copy()
# site = structure[index_in_structure]
new_stru.translate_sites(
indices=np.arange(0, len(structure)), vector=-defect_coord, frac_coords=False, to_unit_cell=False
)
return new_stru
[docs]
def clean_structure(structure: Structure, defect_coord) -> Structure:
"""
Apply successively:
center_wrt_defect(), stru_0_1_to_minus_05_05(), frac_coords_05_to_minus05()
Useful to match structures after clean_structure()
"""
stru = structure.copy()
stru = center_wrt_defect(stru, defect_coord)
stru = stru_0_1_to_minus_05_05(stru)
stru = frac_coords_05_to_minus05(stru)
return stru
[docs]
def map_two_structures_coords(stru_1: Structure, stru_2: Structure, tol: float = 0.1) -> list[int]:
"""
Returns a mapping between two structures based on
the coordinates, within a given tolerance in Angstrom.
stru_1 should be a subset of the superset structure stru_2
"""
cart_1 = stru_1.cart_coords
cart_2 = stru_2.cart_coords
mapping = []
for i, site_1 in enumerate(stru_1): # subset structure
for j, site_2 in enumerate(stru_2): # superset structure
if max(abs(cart_1[i] - cart_2[j])) < tol:
mapping.append(j)
if len(mapping) != len(stru_1):
print(f"Mapping incomplete: only {len(mapping)}/{len(stru_1)} atoms mapped.")
return mapping
[docs]
def accoustic_sum(Hessian_matrix, atom_label) -> np.ndarray:
sum_hessian = np.zeros((3, 3))
for m in range(len(Hessian_matrix[0])):
sum_hessian += Hessian_matrix[m][atom_label]
return sum_hessian
[docs]
def inverse_participation_ratio(eigenvectors):
"""
Compute equation (10) of https://pubs.acs.org/doi/10.1021/acs.chemmater.3c00537
for a given q-point, eigenvectors shape=(nbands,natoms,3)
"""
# for a given q-point, eigenvectors shape=(nbands,natoms,3)
ipr = []
for iband in range(len(eigenvectors)):
sum_atoms = 0
for iatom in range(len(eigenvectors[iband])):
sum_atoms += np.dot(eigenvectors[iband, iatom], eigenvectors[iband, iatom]) ** 2
ipr.append(1 / sum_atoms)
return np.array(ipr).real
[docs]
def localization_ratio(eigenvectors):
"""
See equation (10) of https://pubs.acs.org/doi/10.1021/acs.chemmater.3c00537
for a given q-point, eigenvectors shape=(nbands,natoms,3)
"""
ipr = inverse_participation_ratio(eigenvectors)
return len(eigenvectors[0]) / ipr
[docs]
def vesta_phonon(
eigenvectors,
in_path,
ibands=None,
scale_vector=20,
width_vector=0.3,
color_vector=[255, 0, 0],
centered=True,
factor_keep_vectors=0.1,
out_path="VESTA_FILES",
) -> None:
"""
Draw the phonons eigenvectors on a vesta file.
Inspired from https://github.com/AdityaRoy-1996/Phonopy_VESTA/tree/master
Args:
eigenvectors: phonons eigenvectors for given q-point with shape (nbands,natom,3)
in_path: path where the initial vesta structure in stored, the structure should be consistent
with the eigenvectors provided.
ibands: list of indices of the bands to include, if not provided, all the bands
scale_vector: scaling factor of the vector modulus
width_vector: vector width
color_vector: color in rgb format
centered: center the vesta structure around [0,0,0]
factor_keep_vectors: draw only the eigenvectors with magnitude > factor_keep_vectors * max(magnitude)
out_path: path where .vesta files with vector are stored
"""
vesta = open(in_path).read()
nbands = len(eigenvectors)
natoms = len(eigenvectors[0])
if ibands is None:
ibands = range(nbands)
if os.path.isdir(out_path):
shutil.rmtree(out_path)
os.mkdir(out_path)
else:
os.mkdir(out_path)
path = out_path
for iband in ibands:
towrite = vesta.split("VECTR")[0]
towrite += "VECTR\n"
magnitudes = []
for iatom in range(natoms):
magnitudes.append(
np.sqrt(
eigenvectors[iband][iatom][0] ** 2
+ eigenvectors[iband][iatom][1] ** 2
+ eigenvectors[iband][iatom][2] ** 2
)
)
for iatom in range(natoms):
if magnitudes[iatom] > factor_keep_vectors * max(np.real(magnitudes)):
towrite += "%5d" % (iatom + 1)
towrite += "%10.5f" % (eigenvectors[iband][iatom][0] * int(scale_vector))
towrite += "%10.5f" % (eigenvectors[iband][iatom][1] * int(scale_vector))
towrite += "%10.5f" % (eigenvectors[iband][iatom][2] * int(scale_vector))
towrite += "\n"
towrite += "%5d" % (iatom + 1) + " 0 0 0 0\n 0 0 0 0 0\n"
towrite += "0 0 0 0 0\n"
towrite += "VECTT\n"
for atom in range(natoms):
towrite += "%5d" % (atom + 1)
towrite += f" {width_vector} {color_vector[0]} {color_vector[1]} {color_vector[2]} 0\n"
towrite += "0 0 0 0 0\n"
towrite += "SPLAN"
towrite += vesta.split("SPLAN")[1]
towrite += "VECTS 1.00000"
filename = path + f"/{iband:05}_"
filename += ".vesta"
open(filename, "w").write(towrite)
if centered:
with open(filename) as file:
file_contents = file.read()
search_word = "BOUND\n 0 1 0 1 0 1\n 0 0 0 0 0"
replace_word = (
"BOUND\n -0.5 0.5 -0.5 0.5 -0.5 0.5\n 0 0 0 0 0"
)
updated_contents = file_contents.replace(search_word, replace_word)
with open(filename, "w") as file:
file.write(updated_contents)
print(f"Vesta files created and stored in: \n {os.getcwd()}/{out_path}")