Skip to content

Commit

Permalink
implemented add_dof + necessary symmetry logic
Browse files Browse the repository at this point in the history
  • Loading branch information
@wolearyc
wolearyc committed Jul 29, 2024
1 parent 924cd89 commit 3b6ac3d
Show file tree
Hide file tree
Showing 7 changed files with 372 additions and 278 deletions.
1 change: 1 addition & 0 deletions cspell.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
"dictionaryDefinitions": [],
"dictionaries": [],
"words": [
"argmax",
"AVOGADROS",
"dtype",
"fft",
Expand Down
2 changes: 1 addition & 1 deletion ramannoodle/io/vasp/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
from numpy.typing import NDArray

from ...dynamics import Phonons
from ...polarizability import StructuralSymmetry
from ...symmetry import StructuralSymmetry
from ...globals import ATOMIC_WEIGHTS, ATOMIC_NUMBERS
from .vasp_utils import (
_get_atomic_symbol_from_potcar_line,
Expand Down
222 changes: 71 additions & 151 deletions ramannoodle/polarizability/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,121 +1,18 @@
"""Parent class for polarizability models."""
"""Polarizability models."""

from abc import ABC, abstractmethod

import numpy as np
from numpy.typing import NDArray
from scipy.interpolate import make_interp_spline, BSpline
import spglib

from . import polarizability_utils
from ..symmetry.symmetry_utils import is_orthogonal_to_all
from ..symmetry import StructuralSymmetry
from ..exceptions import InvalidDOFException, SymmetryException


class StructuralSymmetry:
"""Represents symmetries of a crystal structure."""

def __init__( # pylint: disable=too-many-arguments
self,
atomic_numbers: NDArray[np.int32],
lattice: NDArray[np.float64],
fractional_positions: NDArray[np.float64],
symprec: float = 1e-5,
angle_tolerance: float = -1.0,
) -> None:
self._atomic_numbers = atomic_numbers
self._lattice = lattice
self._fractional_positions = fractional_positions

cell = (list(lattice), list(fractional_positions), list(atomic_numbers))
self._symmetry_dict: dict[str, NDArray[np.float64]] | None = (
spglib.get_symmetry(cell, symprec=symprec, angle_tolerance=angle_tolerance)
)
if self._symmetry_dict is None:
raise SymmetryException("symmetry search failed")

self._rotations = self._symmetry_dict["rotations"]
self._translations = self._symmetry_dict["translations"]
self._permutation_matrices = polarizability_utils.compute_permutation_matrices(
self._rotations, self._translations, self._fractional_positions
)

def get_num_nonequivalent_atoms(self) -> int:
"""Returns the number of nonequivalent atoms."""
assert self._symmetry_dict is not None
return len(set(self._symmetry_dict["equivalent_atoms"]))

def get_equivalent_displacements(
self, displacement: NDArray[np.float64]
) -> list[dict[str, list[NDArray[np.float64]]]]:
"""Calculates and returns all symmetrically equivalent displacements.
The return is a little complicated.
[('displacements' : [...], {'transformations' : [...]}]
We want to guarantee that all displacements are orthogonal between
the dictionaries and are all collinear and unique within the
dictionaries.
"""

ref_positions = polarizability_utils.add_fractional_positions(
self._fractional_positions, displacement
)

result = []
orthogonal_displacements: list[NDArray[np.float64]] = []
for rotation, translation, permutation_matrix in zip(
self._rotations, self._translations, self._permutation_matrices
):

# Transform, permute, then get displacements
candidate_positions = polarizability_utils.transform_fractional_positions(
ref_positions, rotation, translation
)
candidate_positions = permutation_matrix @ candidate_positions
candidate_displacement = polarizability_utils.subtract_fractional_positions(
candidate_positions, self._fractional_positions
)

orthogonal_result = polarizability_utils.is_orthogonal_to_all(
candidate_displacement.flatten(),
[item.flatten() for item in orthogonal_displacements],
)

# If new orthogonal displacement discovered
if orthogonal_result == -1:
result.append(
{
"displacements": [candidate_displacement],
"transformations": [(rotation, translation)],
}
)
orthogonal_displacements.append(candidate_displacement)
else:
# Candidate can be collinear to maximum of one orthogonal vector
collinear_result = polarizability_utils.is_non_collinear_with_all(
candidate_displacement,
orthogonal_displacements,
)
if collinear_result != -1: # collinear
collinear_displacements = result[collinear_result]["displacements"]

# Check if we have a duplicate
is_duplicate = False
for collinear_displacement in collinear_displacements:
if np.isclose(
candidate_displacement, collinear_displacement, atol=1e-05
).all():
is_duplicate = True
break
# If unique, add
if not is_duplicate:
result[collinear_result]["displacements"].append(
candidate_displacement
)
result[collinear_result]["transformations"].append(
(rotation, translation)
)

return result


class PolarizabilityModel(ABC): # pylint: disable=too-few-public-methods
"""Represents a polarizability model"""

Expand All @@ -137,12 +34,10 @@ class InterpolationPolarizabilityModel(PolarizabilityModel):
"""

def __init__(
self, structural_symmetry: StructuralSymmetry
) -> None: # - Generate all symmetry equivalent displacements.
def __init__(self, structural_symmetry: StructuralSymmetry) -> None:
self._structural_symmetry = structural_symmetry
self._dof_displacements: list[NDArray[np.float64]] = []
self._dof_interpolations: list[BSpline] = []
self._basis_vectors: list[NDArray[np.float64]] = []
self._interpolations: list[BSpline] = []

def get_polarizability(
self, displacement: NDArray[np.float64]
Expand All @@ -154,55 +49,80 @@ def get_polarizability(

return np.array([])

def add_dof(
def add_dof( # pylint: disable=too-many-locals
self,
displacements_and_polarizabilities: list[
tuple[NDArray[np.float64], NDArray[np.float64]]
],
displacement: NDArray[np.float64],
magnitudes: NDArray[np.float64],
polarizabilities: NDArray[np.float64],
interpolation_dim: int,
) -> None:
"""Adds a degree of freedom as well as any symmetrically
equivalent degrees of freedom."""

if len(displacements_and_polarizabilities) == 0:
raise InvalidDOFException("no dof provided")
"""Adds a degree of freedom"""

dof_displacement, _ = displacements_and_polarizabilities[0]
dof_displacement /= np.linalg.norm(dof_displacement)
parent_displacement = displacement / (np.linalg.norm(displacement) * 10)

# Check that all provided displacements are collinear
for i, (displacement, _) in enumerate(displacements_and_polarizabilities):
if not polarizability_utils.are_collinear(dof_displacement, displacement):
raise InvalidDOFException(f"displacement (index={i}) is not collinear")

# Check that new displacement is orthogonal with existing displacements
result = polarizability_utils.is_orthogonal_to_all(
dof_displacement, self._dof_displacements
)
# Check that the parent displacement is orthogonal to existing basis vectors
result = is_orthogonal_to_all(parent_displacement, self._basis_vectors)
if result != -1:
raise InvalidDOFException(
raise ValueError(
f"new dof is not orthogonal with existing dof (index={result})"
)
if len(magnitudes) == 0:
raise ValueError("no magnitudes provided")
if len(magnitudes) != len(polarizabilities):
raise ValueError(
f"unequal numbers of magnitudes ({len(magnitudes)})) and "
f"polarizabilities ({len(polarizabilities)})"
)

# Get symmetrically equivalent displacements/transformations
self._structural_symmetry.get_equivalent_displacements(dof_displacement)

# Extract any collinear displacements it finds

# Double check that remaining displacements are orthogonal

interpolation = make_interp_spline(
x=[d for d, _ in displacements_and_polarizabilities],
y=[p for _, p in displacements_and_polarizabilities],
k=interpolation_dim,
bc_type="natural",
displacements_and_transformations = (
self._structural_symmetry.get_equivalent_displacements(parent_displacement)
)

# new_dof_displacements, new_dof_interpolations =

self._dof_displacements.append(dof_displacement)
self._dof_interpolations.append(interpolation)
basis_vectors_to_add: list[NDArray[np.float64]] = []
interpolations_to_add: list[BSpline] = []
for dof_dictionary in displacements_and_transformations:
child_displacement = dof_dictionary["displacements"][0]

interpolation_x = []
interpolation_y = []
for collinear_displacement, transformation in zip(
dof_dictionary["displacements"], dof_dictionary["transformations"]
):
_index = np.unravel_index(
np.argmax(np.abs(child_displacement)), child_displacement.shape
)
multiplier = child_displacement[_index] / collinear_displacement[_index]

for magnitude, polarizability in zip(magnitudes, polarizabilities):
interpolation_x.append(multiplier * magnitude)
interpolation_y.append(polarizability @ transformation[0])

# If duplicate magnitudes are generated, too much data has
# been provided
duplicate = polarizability_utils.find_duplicates(interpolation_x)
if duplicate is not None:
raise ValueError(
f"due to symmetry, magnitude {duplicate} should not be specified"
)

if len(interpolation_x) <= interpolation_dim:
raise ValueError(
f"insufficient data ({len(interpolation_x)}) available for"
f"{interpolation_dim}-dimensional interpolation"
)
assert len(interpolation_x) > 1
basis_vectors_to_add.append(
child_displacement / np.linalg.norm(child_displacement)
)
sort_indices = np.argsort(interpolation_x)
interpolations_to_add.append(
make_interp_spline(
x=np.array(interpolation_x)[sort_indices],
y=np.array(interpolation_y)[sort_indices],
k=interpolation_dim,
bc_type=None,
)
)

# # Go through all transformations and fill in for all equivalent atoms
#
self._basis_vectors += basis_vectors_to_add
self._interpolations += interpolations_to_add
Loading

0 comments on commit 3b6ac3d

Please sign in to comment.