reduced third-order invariant of BOP

bop.py

# Copyright (c) Pymatgen Development Team.
from math import acos, sqrt

import numpy as np
from scipy.special import sph_harm
from sympy.physics.wigner import wigner_3j
from pymatgen.analysis.local_env import CrystalNN


def get_neighbor_sites(structure, n, cutoff):
    centsite = structure[n]
    neighsitestmp = [i[0] for i in structure.get_sites_in_sphere(centsite.coords, cutoff)]
    neighsites = []

    if centsite not in neighsitestmp:
        raise ValueError("Could not find center site!")

    neighsitestmp.remove(centsite)
    neighsites = list(neighsitestmp)
    return neighsites


def get_angles(structure, n, neighsites):
    centsite = structure[n]
    centvec = centsite.coords

    thetas = []
    phis = []
    left_of_unity = 1.0 - 1.0e-12

    for j, neigh in enumerate(neighsites):
        rij = neigh.coords - centvec
        vec = rij / np.linalg.norm(rij)
        # z is North pole --> theta between vec and (0, 0, 1)^T.
        # Because vec is normalized, dot product is simply vec[2].
        thetas.append(acos(max(-1.0, min(vec[2], 1.0))))
        tmpphi = 0.0

        # Compute phi only if it is not (almost) perfectly
        # aligned with z-axis.
        if -left_of_unity < vec[2] < left_of_unity:
            # x is prime meridian --> phi between projection of vec
            # into x-y plane and (1, 0, 0)^T
            tmpphi = acos(max(-1.0, min(vec[0] / (sqrt(
                vec[0] * vec[0] + vec[1] * vec[1])), 1.0)))
            if vec[1] < 0.0:
                tmpphi = -tmpphi
        phis.append(tmpphi)

    return thetas, phis


def get_bond_order_parameter(thetas, phis, l=4):
    # (2l+1, len(thetas))
    # Be careful conventions for theta and phi.
    coeffs_lm = np.array([[sph_harm(m, l, phi, theta) for theta, phi in zip(thetas, phis)]
                          for m in range(-l, l + 1)])
    arr_Q_lm = np.average(coeffs_lm, axis=1)
    Q_l_tmp = np.real(np.dot(arr_Q_lm, np.conj(arr_Q_lm)))
    Q_l = np.sqrt(4 * np.pi / (2 * l + 1) * Q_l_tmp)
    return Q_l, arr_Q_lm


def get_third_order_invariant(structure, n, l=4, cutoff=None):
    if cutoff is None:
        cnn = CrystalNN(distance_cutoffs=None, x_diff_weight=0, porous_adjustment=False)
        neighsites = cnn.get_nn(structure, n)
    else:
        neighsites = get_neighbor_sites(structure, n, cutoff)

    thetas, phis = get_angles(structure, n, neighsites)
    if (thetas is None) or (phis is None):
        return None
    Q_l, arr_Q_lm = get_bond_order_parameter(thetas, phis, l)

    W_l = 0

    for m_1 in range(-l, l + 1):
        for m_2 in range(-l, l + 1):
            m_3 = -m_1 - m_2
            if m_3 < -l or m_3 > l:
                continue
            W_l += wigner_3j(l, l, l, m_1, m_2, m_3) \
                * arr_Q_lm[l + m_1] * arr_Q_lm[l + m_2] * arr_Q_lm[l + m_3]

    Q_l_norm_tmp = np.real(np.dot(arr_Q_lm, np.conj(arr_Q_lm)))
    # We need to cast to built-in type
    normalized_W_l = np.real(complex(W_l / (Q_l_norm_tmp ** (3 / 2))))

    return normalized_W_l


if __name__ == '__main__':
    from pymatgen.core import Structure, Lattice
    length_a = 1.0

    fcc_lattice = length_a / np.sqrt(2) * np.array([
        [0, 1, 1],
        [1, 0, 1],
        [1, 1, 0]
    ])
    fcc = Structure(Lattice(fcc_lattice), ['O2-'], [[0, 0, 0]])

    l = 4
    # reduced_W_l = get_third_order_invariant(fcc, 0, l=4, cutoff=length_a)
    reduced_W_l = get_third_order_invariant(fcc, 0, l=4, cutoff=None)
    print(reduced_W_l)