# Source code for hiphive.core.eigentensors

import numpy as np
import itertools
import sympy
from .tensors import rotate_tensor
from ..io.logging import logger
from .utilities import SparseMatrix

logger = logger.getChild('eigensymmetries')

[docs]def create_eigentensors(cs):

for orbit in cs.orbits:
prototype_cluster = cs._cluster_list[orbit.prototype_index]

# The first round of eigentensors are created from the atom label
# permutation condition
ets = init_ets_from_label_symmetry(prototype_cluster)

# Loop over the symmetries and apply them iteratively to reduce the
# number of eigentensors
for rotation_index, permutation_index in orbit.eigensymmetries:

# Rotation amtrix in scaled coordinates
R = cs.rotation_matrices[rotation_index]
p = cs._permutations[permutation_index]  # permutation vector

# Skip if the symmetry is the trivial identity operation
if np.allclose(R, np.eye(3)) and p == tuple(range(orbit.order)):
continue

p_inv = np.argsort(p)

# Create the intermediate constraint matrix and populate it
M = SparseMatrix.zeros(3**orbit.order, len(ets))
populate_constraint_matrix(M, ets, R, p_inv)

nullspace = M.nullspace()

new_ets = []
for solution in nullspace:
# Normalize the solution vector so it is inly integers
# Done by multiplying with the greatest common multiple of the
# denominators present
solution = renormalize_to_integer(solution)
# From the solution vector and the old eigentensors a new et is
# created
new_ets.append(assemble_new_eigentensor(ets, solution))
ets = new_ets

# If the force constant is forbidden by symmetry, break early
if len(ets) == 0:
break

orbit.eigentensors = ets

[docs]def populate_constraint_matrix(M, ets, R, p_inv):
for et_index, et in enumerate(ets):
tmp = rotate_tensor(et.transpose(p_inv), R)
tmp = tmp - et
for element_index, v in enumerate(tmp.flat):
if v != 0:
M[element_index, et_index] = v

[docs]def init_ets_from_label_symmetry(cluster):
""" This function creates eigentensor which fulfill the atom label
symmetry. The only information needed is the corresponding cluster
"""
order = len(cluster)
assert sorted(cluster) == list(cluster)
multiplicities = np.unique(cluster, return_counts=True)[1]
ets = {}
# Loop over all elements, represented by multi indices
for multi_index in itertools.product([0, 1, 2], repeat=order):
# Sort the multi index based on the multiplicities
# e.g. cluster [1 1 2]
# -> [z y x] == [y z x]
sorted_multi_index = []
j = 0
for m in multiplicities:
sorted_multi_index.extend(sorted(multi_index[j:j + m]))
j += m
sorted_multi_index = tuple(sorted_multi_index)
if sorted_multi_index not in ets:
ets[sorted_multi_index] = np.zeros([3]*order, dtype=np.int64)
ets[sorted_multi_index][multi_index] = 1
return list(ets.values())

[docs]def assemble_new_eigentensor(eigentensors, solution):
new_eigentensor = np.zeros(eigentensors[0].shape, dtype=np.int64)
for parameter, eigentensor in zip(solution, eigentensors):
new_eigentensor += np.int64(parameter) * eigentensor
return new_eigentensor

[docs]def renormalize_to_integer(vector):
gcm = 1  # greates common multiple
for v in vector:
if v.is_Integer:
continue
den = v.q  # Denomonator
gcd = sympy.gcd(gcm, den)  # greates common divisor
tmp = den / gcd
gcm = gcm * tmp
# Now all elements should be smallest possible integers
vector = vector * gcm
return vector