# OctaDist Copyright (C) 2019-2024 Rangsiman Ketkaew et al.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# 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 tkinter as tk
from tkinter import scrolledtext as tkscrolled
import numpy as np
import rmsd
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from octadist.src import linear, elements, plane, util
[docs]
class CalcJahnTeller:
"""
Calculate angular Jahn-Teller distortion parameter [1]_.
Parameters
----------
atom : array_like
Atomic labels of full complex.
coord : array_like
Atomic coordinates of full complex.
master : None, object
If None, use tk.Tk().
If not None, use tk.Toplevel(master).
cutoff_global : int or float
Global cutoff for screening bonds.
Default is 2.0.
cutoff_hydrogen : int or float
Cutoff for screening hydrogen bonds.
Default is 1.2.
icon : str, optional
If None, use tkinter default icon.
If not None, use user-defined icon.
Examples
--------
>>> atom = ['Fe', 'N', 'N', 'N', 'O', 'O', 'O']
>>> coord = [[2.298354000, 5.161785000, 7.971898000],
[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 = CalcJahnTeller(atom=atom, coord=coord)
>>> test.start_app()
>>> test.find_bond()
>>> test.show_app()
References
----------
.. [1] J. M. Holland, J. A. McAllister, C. A. Kilner,
M. Thornton-Pett, A. J. Bridgeman, M. A. Halcrow.
Stereochemical effects on the spin-state transition
shown by salts of [FeL2]2+ [L = 2,6-di(pyrazol-1-yl)pyridine].
J. Chem. Soc., Dalton Trans., 2002, 548-554.
DOI: 10.1039/B108468M.
"""
def __init__(
self,
atom,
coord,
cutoff_global=2.0,
cutoff_hydrogen=1.2,
master=None,
icon=None,
):
self.atom = atom
self.coord = coord
self.cutoff_global = cutoff_global
self.cutoff_hydrogen = cutoff_hydrogen
if master is None:
self.wd = tk.Tk()
else:
self.wd = tk.Toplevel(master)
self.icon = icon
self.bond_list = []
self.coord_A = []
self.coord_B = []
[docs]
def start_app(self):
"""
Start application.
"""
if self.icon is not None:
self.wd.wm_iconbitmap(self.icon)
self.wd.title("Calculate Jahn-Teller distortion parameter")
self.wd.geometry("630x550")
self.wd.resizable(0, 0)
self.lbl = tk.Label(self.wd, text="Group A")
self.lbl.config(width=12)
self.lbl.grid(padx="10", pady="5", row=0, column=0, columnspan=2)
self.lbl = tk.Label(self.wd, text="Group B")
self.lbl.config(width=12)
self.lbl.grid(padx="10", pady="5", row=0, column=2, columnspan=2)
self.box_1 = tkscrolled.ScrolledText(
self.wd, height="12", width="40", wrap="word", undo="True"
)
self.box_1.grid(padx="5", pady="5", row=1, column=0, columnspan=2)
self.box_2 = tkscrolled.ScrolledText(
self.wd, height="12", width="40", wrap="word", undo="True"
)
self.box_2.grid(padx="5", pady="5", row=1, column=2, columnspan=2)
self.btn = tk.Button(
self.wd,
text="Select ligand set A",
command=lambda: self.pick_atom(group="A"),
)
self.btn.config(width=15, relief=tk.RAISED)
self.btn.grid(padx="10", pady="5", row=2, column=0, columnspan=2)
self.btn = tk.Button(
self.wd,
text="Select ligand set B",
command=lambda: self.pick_atom(group="B"),
)
self.btn.config(width=15, relief=tk.RAISED)
self.btn.grid(padx="10", pady="5", row=2, column=2, columnspan=2)
self.btn = tk.Button(
self.wd, text="Calculate parameter", command=lambda: self.plot_fit_plane()
)
self.btn.config(width=15, relief=tk.RAISED)
self.btn.grid(padx="10", pady="5", row=3, column=0, columnspan=2)
self.btn = tk.Button(
self.wd, text="Clear all", command=lambda: self.clear_text()
)
self.btn.config(width=15, relief=tk.RAISED)
self.btn.grid(padx="10", pady="5", row=3, column=2, columnspan=2)
self.lbl = tk.Label(
self.wd, text="Supplementary angles between two planes (in degree)"
)
self.lbl.grid(pady="10", row=4, columnspan=4)
self.lbl_angle1 = tk.Label(self.wd, text="Angle 1")
self.lbl_angle1.grid(pady="5", row=5, column=0)
self.box_angle1 = tk.Entry(self.wd, width="20", justify="center")
self.box_angle1.grid(row=5, column=1, sticky=tk.W)
self.lbl_angle2 = tk.Label(self.wd, text="Angle 2")
self.lbl_angle2.grid(pady="5", row=6, column=0)
self.box_angle2 = tk.Entry(self.wd, width="20", justify="center")
self.box_angle2.grid(row=6, column=1, sticky=tk.W)
self.lbl = tk.Label(self.wd, text="The equation of the planes")
self.lbl.grid(pady="10", row=7, columnspan=4)
self.lbl_eq1 = tk.Label(self.wd, text="Plane A ")
self.lbl_eq1.grid(pady="5", row=8, column=0)
self.box_eq1 = tk.Entry(self.wd, width="60", justify="center")
self.box_eq1.grid(pady="5", row=8, column=1, columnspan=2, sticky=tk.W)
self.lbl_eq2 = tk.Label(self.wd, text="Plane B ")
self.lbl_eq2.grid(pady="5", row=9, column=0)
self.box_eq2 = tk.Entry(self.wd, width="60", justify="center")
self.box_eq2.grid(pady="5", row=9, column=1, columnspan=2, sticky=tk.W)
[docs]
def find_bond(self):
"""
Find bonds.
See Also
--------
octadist.src.util.find_bonds :
Find atomic bonds.
"""
_, self.bond_list = util.find_bonds(
self.atom, self.coord, self.cutoff_global, self.cutoff_hydrogen
)
#################
# Picking atoms #
#################
[docs]
def pick_atom(self, group):
"""
On-mouse pick atom and get XYZ coordinate.
Parameters
----------
group : str
Group A or B.
"""
fig = plt.figure()
# fig = plt.figure(figsize=(5, 4), dpi=100)
ax = Axes3D(fig)
# ax = fig.add_subplot(111, projection='3d')
# Plot all atoms
for i in range(len(self.coord)):
# Determine atomic number
n = elements.number_to_symbol(self.atom[i])
ax.scatter(
self.coord[i][0],
self.coord[i][1],
self.coord[i][2],
marker="o",
linewidths=0.5,
edgecolors="black",
picker=5,
color=elements.number_to_color(n),
label=f"{self.atom[i]}",
s=elements.number_to_radii(n) * 300,
)
atoms_pair = []
for i in range(len(self.bond_list)):
get_atoms = self.bond_list[i]
x, y, z = zip(*get_atoms)
atoms = list(zip(x, y, z))
atoms_pair.append(atoms)
# Draw line
for i in range(len(atoms_pair)):
merge = list(zip(atoms_pair[i][0], atoms_pair[i][1]))
x, y, z = merge
ax.plot(x, y, z, "k-", color="black", linewidth=2)
# Set legend
# Remove duplicate labels in legend.
# Ref.https://stackoverflow.com/a/26550501/6596684
handles, labels = ax.get_legend_handles_labels()
handle_list, label_list = [], []
for handle, label in zip(handles, labels):
if label not in label_list:
handle_list.append(handle)
label_list.append(label)
leg = fig.legend(
handle_list, label_list, loc="lower left", scatterpoints=1, fontsize=12
)
# Fixed size of point in legend
# Ref. https://stackoverflow.com/a/24707567/6596684
for i in range(len(leg.legendHandles)):
leg.legendHandles[i]._sizes = [90]
# Set axis
ax.set_xlabel(r"X", fontsize=15)
ax.set_ylabel(r"Y", fontsize=15)
ax.set_zlabel(r"Z", fontsize=15)
ax.set_title("Full complex", fontsize="12")
ax.grid(True)
def insert_text(text, coord, group):
"""
Insert text in boxes.
Parameters
----------
text : str
Text.
coord : list or array
Coordinates.
group : str
Group A or B.
"""
if group == "A":
self.box_1.insert(tk.INSERT, text + "\n")
self.box_1.see(tk.END)
self.coord_A.append(coord)
elif group == "B":
self.box_2.insert(tk.INSERT, text + "\n")
self.box_2.see(tk.END)
self.coord_B.append(coord)
def on_pick(event):
"""
Pick point and get XYZ data
Parameters
----------
event : object
Event object for on-pick function.
Examples
--------
>>> def on_pick(event):
>>> ... ind = event.ind
>>> ... print("on_pick scatter:", ind, np.take(x, ind), np.take(y, ind))
"""
ind = event.ind[0]
x, y, z = event.artist._offsets3d
for i in range(len(self.atom)):
if x[ind] == self.coord[i][0]:
if y[ind] == self.coord[i][1]:
if z[ind] == self.coord[i][2]:
atom = self.atom[i]
break
results = f"{i + 1} {atom} : {x[ind]} {y[ind]} {z[ind]}"
coord = [x[ind], y[ind], z[ind]]
insert_text(results, coord, group)
# Highlight selected atom
index = elements.number_to_symbol(atom)
ax.scatter(
x[ind],
y[ind],
z[ind],
marker="o",
linewidths=0.5,
edgecolors="orange",
picker=5,
alpha=0.5,
color="yellow",
s=elements.number_to_radii(index) * 400,
)
# print(i+1, atom, x[ind], y[ind], z[ind])
fig.canvas.mpl_connect("pick_event", on_pick)
# plt.axis('equal')
# plt.axis('off')
plt.show()
########################################
# Plot fit plant to the selected atoms #
########################################
[docs]
def plot_fit_plane(self):
"""
Display complex and two fit planes of two sets of ligand in molecule.
"""
###############
# Clear boxes #
###############
self.box_angle1.delete(0, tk.END)
self.box_angle2.delete(0, tk.END)
self.box_eq1.delete(0, tk.END)
self.box_eq2.delete(0, tk.END)
########################
# Find eq of the plane #
########################
xx, yy, z, abcd = plane.find_fit_plane(self.coord_A)
plane_a = (xx, yy, z)
a1, b1, c1, d1 = abcd
self.box_eq1.insert(
tk.INSERT, f"{a1:8.5f}x {b1:+8.5f}y {c1:+8.5f}z {d1:+8.5f} = 0"
)
xx, yy, z, abcd = plane.find_fit_plane(self.coord_B)
plane_b = (xx, yy, z)
a2, b2, c2, d2 = abcd
self.box_eq2.insert(
tk.INSERT, f"{a2:8.5f}x {b2:+8.5f}y {c2:+8.5f}z {d2:+8.5f} = 0"
)
####################################
# Calculate angle between 2 planes #
####################################
angle = linear.angle_btw_planes(a1, b1, c1, a2, b2, c2)
self.box_angle1.insert(tk.INSERT, f"{angle:10.6f}") # insert to box
sup_angle = abs(180 - angle) # supplementary angle
self.box_angle2.insert(tk.INSERT, f"{sup_angle:10.6f}") # insert to box
###############
# Plot planes #
###############
fig = plt.figure()
# fig = plt.figure(figsize=(5, 4), dpi=100)
ax = Axes3D(fig)
# ax = fig.add_subplot(111, projection='3d')
# Plot all atoms
for i in range(len(self.coord)):
# Determine atomic number
n = elements.number_to_symbol(self.atom[i])
ax.scatter(
self.coord[i][0],
self.coord[i][1],
self.coord[i][2],
marker="o",
linewidths=0.5,
edgecolors="black",
picker=5,
color=elements.number_to_color(n),
label=f"{self.atom[i]}",
s=elements.number_to_radii(n) * 300,
)
atoms_pair = []
for i in range(len(self.bond_list)):
get_atoms = self.bond_list[i]
x, y, z = zip(*get_atoms)
atoms = list(zip(x, y, z))
atoms_pair.append(atoms)
# Draw line
for i in range(len(atoms_pair)):
merge = list(zip(atoms_pair[i][0], atoms_pair[i][1]))
x, y, z = merge
ax.plot(x, y, z, "k-", color="black", linewidth=2)
# Set legend
# Remove duplicate labels in legend.
# Ref.https://stackoverflow.com/a/26550501/6596684
handles, labels = ax.get_legend_handles_labels()
handle_list, label_list = [], []
for handle, label in zip(handles, labels):
if label not in label_list:
handle_list.append(handle)
label_list.append(label)
leg = fig.legend(
handle_list, label_list, loc="lower left", scatterpoints=1, fontsize=12
)
# Fixed size of point in legend
# Ref. https://stackoverflow.com/a/24707567/6596684
for i in range(len(leg.legendHandles)):
leg.legendHandles[i]._sizes = [90]
# Set axis
ax.set_xlabel(r"X", fontsize=15)
ax.set_ylabel(r"Y", fontsize=15)
ax.set_zlabel(r"Z", fontsize=15)
ax.set_title("Full complex", fontsize="12")
ax.grid(True)
# Plot plane A
xx, yy, z = plane_a
ax.plot_surface(xx, yy, z, alpha=0.2, color="green")
# Plot plane B
xx, yy, z = plane_b
ax.plot_surface(xx, yy, z, alpha=0.2, color="red")
# ax.set_xlim(-10, 10)
# ax.set_ylim(-10, 10)
# ax.set_zlim(0, 10)
# plt.axis('equal')
# plt.axis('off')
plt.show()
##################
# Clear text box #
##################
[docs]
def clear_text(self):
"""
Clear text in box A & B.
"""
self.box_1.delete(1.0, tk.END)
self.box_2.delete(1.0, tk.END)
self.box_angle1.delete(0, tk.END)
self.box_angle2.delete(0, tk.END)
self.box_eq1.delete(0, tk.END)
self.box_eq2.delete(0, tk.END)
self.coord_A = []
self.coord_B = []
try:
plt.close("all")
except AttributeError:
pass
[docs]
def show_app(self):
"""
Show application.
"""
self.wd.mainloop()
[docs]
class CalcRMSD:
"""
Calculate root mean squared displacement of atoms in complex, RMSD [2]_.
Parameters
----------
coord_1 : array_like
Atomic coordinates of structure 1.
coord_2 : array_like
Atomic coordinates of structure 2.
atom_1 : list or tuple, optional
Atomic symbols of structure 1.
atom_2 : list or tuple, optional
Atomic symbols of structure 2.
If no atom_2 specified, assign it with None.
Returns
-------
rmsd_normal : float
Normal RMSD.
rmsd_translate : float
Translate RMSD (re-centered).
rmsd_rotate : float
Kabsch RMSD (rotated).
References
----------
.. [2] J. C. Kromann. https://github.com/charnley/rmsd.
Examples
--------
>>> # Example of structure 1
>>> comp1 = [[10.1873, 5.7463, 5.615],
[8.494, 5.9735, 4.8091],
[9.6526, 6.4229, 7.3079],
[10.8038, 7.5319, 5.1762],
[9.6229, 3.9221, 6.0083],
[12.0065, 5.5562, 6.3497],
[10.8046, 4.9471, 3.9219]]
>>> # Example of structure 1
>>> comp2 = [[12.0937, 2.4505, 3.4207],
[12.9603, 2.2952, 1.7286],
[13.4876, 1.6182, 4.4230],
[12.8522, 4.3174, 3.9894],
[10.9307, 0.7697, 2.9315],
[10.7878, 2.2987, 5.1071],
[10.6773, 3.7960, 2.5424]]
>>> test = CalcRMSD(coord_1=comp1, coord_2=comp2)
>>> test.calc_rmsd()
>>> test.rmsd_normal
6.758144
>>> test.rmsd_translate
0.305792
>>> test.rmsd_rotate
0.277988
"""
def __init__(
self, coord_1, coord_2, atom_1=None, atom_2=None, master=None, icon=None
):
self.coord_1 = np.asarray(coord_1, dtype=np.float64)
self.coord_2 = np.asarray(coord_2, dtype=np.float64)
self.atom_1 = atom_1
self.atom_2 = atom_2
if self.atom_2 is None:
self.atom_2 = self.atom_1
if master is None:
self.wd = tk.Tk()
else:
self.wd = tk.Toplevel(master)
self.icon = icon
self.rmsd_normal = 0
self.rmsd_translate = 0
self.rmsd_rotate = 0
[docs]
def start_app(self):
if self.icon is not None:
self.wd.wm_iconbitmap(self.icon)
self.wd.title("Calculate RMSD")
self.wd.geometry("700x440")
self.wd.resizable(0, 0)
self.lbl = tk.Label(
self.wd, text="Calculate root-mean-square deviation of atomic positions"
)
# self.lbl.config(width=12)
self.lbl.grid(padx="5", pady="5", row=0, column=0, columnspan=2)
self.lbl = tk.Label(self.wd, text="Structure A")
self.lbl.config(width=12)
self.lbl.grid(padx="10", pady="5", row=1, column=0)
self.lbl = tk.Label(self.wd, text="Structure B")
self.lbl.config(width=12)
self.lbl.grid(padx="10", pady="5", row=1, column=1)
self.box_1 = tkscrolled.ScrolledText(
self.wd, height="12", width="45", wrap="word", undo="True"
)
self.box_1.grid(padx="5", pady="5", row=2, column=0)
self.box_2 = tkscrolled.ScrolledText(
self.wd, height="12", width="45", wrap="word", undo="True"
)
self.box_2.grid(padx="5", pady="5", row=2, column=1)
self.lbl = tk.Label(self.wd, text="Normal RMSD")
self.lbl.grid(padx="5", pady="5", row=3, column=0, sticky=tk.E)
self.box_3 = tk.Entry(self.wd, width="20", justify="center")
self.box_3.grid(padx="5", row=3, column=1, sticky=tk.W)
self.lbl = tk.Label(self.wd, text="Re-centered RMSD")
self.lbl.grid(padx="5", pady="5", row=4, column=0, sticky=tk.E)
self.box_4 = tk.Entry(self.wd, width="20", justify="center")
self.box_4.grid(padx="5", row=4, column=1, sticky=tk.W)
self.lbl = tk.Label(self.wd, text="Rotated RMSD")
self.lbl.grid(padx="5", pady="5", row=5, column=0, sticky=tk.E)
self.box_5 = tk.Entry(self.wd, width="20", justify="center")
self.box_5.grid(padx="5", row=5, column=1, sticky=tk.W)
self.btn = tk.Button(
self.wd, text="Calculate RMSD", command=lambda: self.calc_and_show()
)
self.btn.config(width=15, relief=tk.RAISED)
self.btn.grid(padx="10", pady="15", row=6, column=0, columnspan=2)
# When app starts, show atomic coordinates in box.
self.show_coord()
[docs]
def show_coord(self):
"""
Show atomic coordinates in box.
"""
num_atoms = len(self.coord_1)
self.box_1.insert(tk.INSERT, f"Total atoms: {num_atoms}\n\n")
self.box_2.insert(tk.INSERT, f"Total atoms: {num_atoms}\n\n")
if self.atom_1 is None or self.atom_2 is None:
self.atom_1 = ["N/A"] * num_atoms
self.atom_2 = ["N/A"] * num_atoms
for i in range(num_atoms):
self.box_1.insert(
tk.INSERT, f"{i+1}. {self.atom_1[i]}\t{self.coord_1[i]}\n"
)
self.box_2.insert(
tk.INSERT, f"{i+1}. {self.atom_2[i]}\t{self.coord_2[i]}\n"
)
[docs]
def calc_rmsd(self):
"""
Calculate normal, translated, and rotated RMSD.
"""
tmp_1, tmp_2 = self.coord_1.copy(), self.coord_2.copy()
self.rmsd_normal = rmsd.rmsd(self.coord_1, self.coord_2)
# Manipulate recenter
self.coord_1 -= rmsd.centroid(self.coord_1)
self.coord_2 -= rmsd.centroid(self.coord_2)
self.rmsd_translate = rmsd.rmsd(self.coord_1, self.coord_2)
# Rotate
rotation_matrix = rmsd.kabsch(self.coord_1, self.coord_2)
self.coord_1 = np.dot(self.coord_1, rotation_matrix)
self.rmsd_rotate = rmsd.rmsd(self.coord_1, self.coord_2)
self.coord_1, self.coord_2 = tmp_1, tmp_2
[docs]
def calc_and_show(self):
"""
Execute calc_rmsd function to calculate RMSD
and show results in box.
"""
self.calc_rmsd()
self.box_3.delete(0, tk.END)
self.box_4.delete(0, tk.END)
self.box_5.delete(0, tk.END)
self.box_3.insert(tk.INSERT, f"{self.rmsd_normal: 3.6f}")
self.box_4.insert(tk.INSERT, f"{self.rmsd_translate: 3.6f}")
self.box_5.insert(tk.INSERT, f"{self.rmsd_rotate: 3.6f}")
[docs]
def show_app(self):
"""
Show application.
"""
self.wd.mainloop()