# OctaDist  Copyright (C) 2019  Rangsiman Ketkaew et al.
#
# This program is free software: you can redistribute it and/or modify
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import numpy as np
from scipy.spatial import distance

from octadist.src import linear, plane, projection

[docs]class CalcDistortion:
"""
Calculate octahedral histortion parameters:

- Bond distance : :meth:calc_d_bond
- Mean bond distance : :meth:calc_d_mean
- Bond angle around metal center atom : :meth:calc_bond_angle
- zeta parameter : :meth:calc_zeta
- Delta parameter : :meth:calc_delta
- Sigma parameter : :meth:calc_sigma
- Minimum Tehta parameter : :meth:calc_theta_min
- Maximum Theta parameter : :meth:calc_theta_max
- Mean Theta parametes : :meth:calc_theta

Parameters
----------
coord : array_like
Atomic coordinates of octahedral structure.

Examples
--------
>>> coord = [[2.298354000, 5.161785000, 7.971898000],  # <- Metal atom
[1.885657000, 4.804777000, 6.183726000],
[1.747515000, 6.960963000, 7.932784000],
[4.094380000, 5.807257000, 7.588689000],
[0.539005000, 4.482809000, 8.460004000],
[2.812425000, 3.266553000, 8.131637000],
[2.886404000, 5.392925000, 9.848966000]]
>>> test = CalcDistortion(coord)
>>> test.sigma
47.926528379270124

"""

def __init__(self, coord):
self.coord = coord

if type(self.coord) == np.ndarray:
pass
else:
self.coord = np.asarray(self.coord, dtype=np.float64)

self.bond_dist = []
self.d_mean = 0
self.diff_dist = []
self.zeta = 0
self.delta = 0
self.cis_angle = 0
self.trans_angle = 0
self.sigma = 0
self.eight_theta = []
self.theta = 0
self.theta_min = 0
self.theta_max = 0
self.eq_of_plane = []
self.non_octa = False

self.calc_d_bond()
self.calc_d_mean()
self.calc_bond_angle()
self.calc_zeta()
self.calc_delta()
self.calc_sigma()
self.calc_theta()
self.calc_theta_min()
self.calc_theta_max()

[docs]    def calc_d_bond(self):
"""
Calculate metal-ligand bond distance and return value in Angstrom.

--------
calc_d_mean :
Calculate mean metal-ligand bond length.

"""
self.bond_dist = [distance.euclidean(self.coord[0], self.coord[i]) for i in range(1, 7)]
self.bond_dist = np.asarray(self.bond_dist, dtype=np.float64)

[docs]    def calc_d_mean(self):
"""
Calculate mean distance parameter and return value in Angstrom.

--------
calc_d_bond :
Calculate metal-ligand bonds length.

"""
self.d_mean = np.mean(self.bond_dist)

[docs]    def calc_bond_angle(self):
"""
Calculate 12 cis and 3 trans unique angles in octahedral structure.

--------
calc_sigma :
Calculate Sigma parameter.

"""
all_angle = []
for i in range(1, 7):
for j in range(i + 1, 7):
vec1 = self.coord[i] - self.coord[0]
vec2 = self.coord[j] - self.coord[0]
angle = linear.angle_btw_vectors(vec1, vec2)
all_angle.append(angle)

# Sort the angle from the lowest to the highest
sorted_angle = sorted(all_angle)
self.cis_angle = [sorted_angle[i] for i in range(12)]
self.trans_angle = [sorted_angle[i] for i in range(12, 15)]

[docs]    def calc_zeta(self):
"""
Calculate zeta parameter [1]_ and return value in Angstrom.

--------
calc_d_bond :
Calculate metal-ligand bonds length.
calc_d_mean :
Calculate mean metal-ligand bond length.

References
----------
.. [1] M. Buron-Le Cointe, J. Hébert, C. Baldé, N. Moisan,
L. Toupet, P. Guionneau, J. F. Létard, E. Freysz,
H. Cailleau, and E. Collet. - Intermolecular control of
thermoswitching and photoswitching phenomena in two
spin-crossover polymorphs. Phys. Rev. B 85, 064114.

"""
diff_dist = [abs(self.bond_dist[i] - self.d_mean) for i in range(6)]
self.diff_dist = np.asarray(diff_dist, dtype=np.float64)

self.zeta = np.sum(self.diff_dist)

[docs]    def calc_delta(self):
"""
Calculate Delta parameter, also known as Tilting distortion parameter [2]_.

--------
calc_d_bond :
Calculate metal-ligand bonds length.
calc_d_mean :
Calculate mean metal-ligand bond length.

References
----------
.. [2] M. W. Lufaso and P. M. Woodward. - Jahn–Teller distortions,
cation ordering and octahedral tilting in perovskites.
Acta Cryst. (2004). B60, 10-20. DOI: 10.1107/S0108768103026661

"""
delta = sum(pow((self.bond_dist[i] - self.d_mean) / self.d_mean, 2) for i in range(6))
self.delta = delta / 6

[docs]    def calc_sigma(self):
"""
Calculate Sigma parameter [3]_ and return value in degree.

--------
calc_bond_angle :
Calculate bond angles between ligand-metal-ligand.

References
----------
.. [3] James K. McCusker, A. L. Rheingold, D. N. Hendrickson.
Variable-Temperature Studies of Laser-Initiated 5T2 → 1A1
Intersystem Crossing in Spin-Crossover Complexes:
Empirical Correlations between Activation Parameters
and Ligand Structure in a Series of Polypyridyl.
Ferrous Complexes. Inorg. Chem. 1996, 35, 2100.

"""
self.sigma = sum(abs(90.0 - self.cis_angle[i]) for i in range(12))

[docs]    def determine_faces(self):
"""
Refine the order of ligand atoms in order to find the plane for projection.

Returns
-------
coord_metal : array_like
Coordinate of metal atom.
coord_lig : array_like
Coordinate of ligand atoms.

--------
calc_theta :
Calculate mean Theta parameter

Examples
--------
>>> bef = np.array([
[4.0674, 7.2040, 13.6117]
[4.3033, 7.3750, 11.7292]
[3.8326, 6.9715, 15.4926]
[5.8822, 6.4461, 13.4312]
[3.3002, 5.3828, 13.6316]
[4.8055, 8.9318, 14.2716]
[2.3184, 8.0165, 13.1152]
])
>>> metal, coord = self.determine_faces(bef)
>>> metal
[ 4.0674  7.204  13.6117]
>>> coord_lig
[[ 4.3033  7.375  11.7292]      # Front face
[ 4.8055  8.9318 14.2716]      # Front face
[ 5.8822  6.4461 13.4312]      # Front face
[ 2.3184  8.0165 13.1152]      # Back face
[ 3.8326  6.9715 15.4926]      # Back face
[ 3.3002  5.3828 13.6316]]     # Back face

"""
# Metal and ligand atoms
coord_metal = self.coord[0]
ligands = self.coord[1:]
coord_lig = np.array([self.coord[i] for i in range(1, 7)])

# Find vector from metal to ligand atoms
metal_to_lig = coord_lig - coord_metal

# Find maximum angle
max_angle = self.trans_angle[0]

# Identify which N is in line with ligand 1
def_change = 6
for n in range(6):
test = linear.angle_btw_vectors(metal_to_lig[0], metal_to_lig[n])
if test > (max_angle - 1):
def_change = n

test_max = 0
new_change = 0
for n in range(6):
test = linear.angle_btw_vectors(metal_to_lig[0], metal_to_lig[n])
if test > test_max:
test_max = test
new_change = n

# Check if the structure is octahedron or not
if def_change != new_change:
self.non_octa = True
def_change = new_change

# Swap ligand
# As ligands in 2-dimensions array, we use the following technique
# to swap two lists of elements with avoiding passing reference value
ligands[[4, def_change]] = ligands[[def_change, 4]]

# Update vector from metal to ligand atoms
metal_to_lig = ligands - coord_metal

# Identify which N is in line with ligand 2
for n in range(6):
test = linear.angle_btw_vectors(metal_to_lig[1], metal_to_lig[n])
if test > (max_angle - 1):
def_change = n

test_max = 0
for n in range(6):
test = linear.angle_btw_vectors(metal_to_lig[1], metal_to_lig[n])
if test > test_max:
test_max = test
new_change = n

if def_change != new_change:
self.non_octa = True
def_change = new_change

# Swap ligand
ligands[[5, def_change]] = ligands[[def_change, 5]]

# Update vector from metal to ligand atoms
metal_to_lig = ligands - coord_metal

# Identify which N is in line with ligand 3
for n in range(6):
test = linear.angle_btw_vectors(metal_to_lig[2], metal_to_lig[n])
if test > (max_angle - 1):
def_change = n

test_max = 0
for n in range(6):
test = linear.angle_btw_vectors(metal_to_lig[2], metal_to_lig[n])
if test > test_max:
test_max = test
new_change = n

if def_change != new_change:
self.non_octa = True
def_change = new_change

# Swap ligand
ligands[[3, def_change]] = ligands[[def_change, 3]]

# New atom order
coord_lig = np.array([ligands[i] for i in range(6)])

return coord_metal, coord_lig

[docs]    def calc_theta(self):
"""
Calculate Theta parameter [4]_ and value in degree.

--------
calc_theta_min :
Calculate minimum Theta parameter.
calc_theta_max :
Calculate maximum Theta parameter.
Calculate cosine angle between two vectors.
Calculate cosine angle between two vectors sensitive to CW/CCW direction.
Find the equation of the plane.
Orthogonal projection of point onto the plane.

References
----------
.. [4] M. Marchivie, P. Guionneau, J.-F. Létard, D. Chasseau.
Photo‐induced spin‐transition: the role of the iron(II)
environment distortion. Acta Crystal-logr. Sect. B Struct.
Sci. 2005, 61, 25.

"""
# Get refined atomic coordinates
coord_metal, coord_lig = self.determine_faces()

# loop over 8 faces
for r in range(8):
a, b, c, d = plane.find_eq_of_plane(coord_lig[0], coord_lig[1], coord_lig[2])
self.eq_of_plane.append([a, b, c, d])

# Project metal and other three ligand atom onto the plane
projected_m = projection.project_atom_onto_plane(coord_metal, a, b, c, d)
projected_lig4 = projection.project_atom_onto_plane(coord_lig[3], a, b, c, d)
projected_lig5 = projection.project_atom_onto_plane(coord_lig[4], a, b, c, d)
projected_lig6 = projection.project_atom_onto_plane(coord_lig[5], a, b, c, d)

# Find the vectors between atoms that are on the same plane
# These vectors will be used to calculate Theta afterward.
vector_theta = np.array(
[
coord_lig[0] - projected_m,
coord_lig[1] - projected_m,
coord_lig[2] - projected_m,
projected_lig4 - projected_m,
projected_lig5 - projected_m,
projected_lig6 - projected_m,
]
)

# Check if the direction is CW or CCW
a12 = linear.angle_btw_vectors(vector_theta[0], vector_theta[1])
a13 = linear.angle_btw_vectors(vector_theta[0], vector_theta[2])

# If angle of interest is smaller than its neighbor,
# define it as CW direction, if not, it will be CCW instead.
if a12 < a13:
direction = np.cross(vector_theta[0], vector_theta[1])
else:
direction = np.cross(vector_theta[2], vector_theta[0])

# Calculate individual theta angle
theta1 = linear.angle_sign(vector_theta[0], vector_theta[3], direction)
theta2 = linear.angle_sign(vector_theta[3], vector_theta[1], direction)
theta3 = linear.angle_sign(vector_theta[1], vector_theta[4], direction)
theta4 = linear.angle_sign(vector_theta[4], vector_theta[2], direction)
theta5 = linear.angle_sign(vector_theta[2], vector_theta[5], direction)
theta6 = linear.angle_sign(vector_theta[5], vector_theta[0], direction)

indi_theta = np.array([theta1, theta2, theta3, theta4, theta5, theta6])

self.eight_theta.append(sum(abs(indi_theta - 60)))

# Use deep copy so as to avoid pass by reference
tmp = coord_lig[1].copy()
coord_lig[1] = coord_lig[3].copy()
coord_lig[3] = coord_lig[5].copy()
coord_lig[5] = coord_lig[2].copy()
coord_lig[2] = tmp.copy()

# If 3rd round, permutation face will be switched
# from N1N2N3
# to N1N4N2,
# to N1N6N4,
# to N1N3N6, and then back to N1N2N3
if r == 3:
coord_lig[[0, 4]] = coord_lig[[4, 0]]
coord_lig[[1, 5]] = coord_lig[[5, 1]]
coord_lig[[2, 3]] = coord_lig[[3, 2]]

self.theta = sum(self.eight_theta) / 2

[docs]    def calc_theta_min(self):
"""
Calculate minimum Theta parameter and return value in degree.

--------
calc_theta :
Calculate mean Theta parameter

"""
sorted_theta = sorted(self.eight_theta)
self.theta_min = sum(sorted_theta[i] for i in range(4))

[docs]    def calc_theta_max(self):
"""
Calculate maximum Theta parameter and return value in degree.