"""
Molecular geometry tools built on top of ASE: alignment, transformation,
PCA-based matching, and VMD-based connectivity extraction.
"""
from ase.calculators.gaussian import Gaussian
import numpy as np
from ase import Atom
from eproducer.miscellaneous import output_terminal
[docs]
class MoleculeSetter:
"""
Tool to edit or extract information from a molecule.
Parameters
==========
atoms: ase.Atoms
molecule to be analyzed.
"""
def __init__(self, atoms):
"""
Stores the molecule to be edited. See the class docstring for the
parameters.
"""
self.atoms = atoms
[docs]
def rot_x(self, angle):
"""
Rotation matrix around x axis.
Parameters
==========
angle: float[radians]
angle to rotate around the x axis
Return
======
(numpy.array) [3float x 3float] Rotation matrix.
"""
c = np.cos(angle)
s = np.sin(angle)
R = np.array([[1, 0, 0],
[0, c, -s],
[0, s, c]])
return R
[docs]
def rot_y(self, angle):
"""
Retuns the rotation matrix around y axis.
Parameters
==========
angle: float[radians]
angle to rotate around the y axis
Return
======
(numpy.array) [3float x 3float] Rotation matrix.
"""
c = np.cos(angle)
s = np.sin(angle)
R = np.array([[c, 0, s],
[0, 1, 0],
[-s, 0, c]])
return R
[docs]
def rot_z(self, angle):
"""
Retuns the rotation matrix around z axis.
Parameters
==========
angle: float[radians]
angle to rotate around the z axis
Return
======
(numpy.array) [3float x 3float] Rotation matrix.
"""
c = np.cos(angle)
s = np.sin(angle)
R = np.array([[c, -s, 0],
[s, c, 0],
[0, 0, 1]])
return R
[docs]
def align_axis(self, vector):
"""
Computes the rotation matrix that aligns a vector with the positive x
axis. The atoms are not modified.
Parameters
==========
vector: array
vector to be aligned.
Return
======
(numpy.array) [3float x 3float] Transformation matrix.
"""
xyproj = vector.copy()
xyproj[2] = 0
phi = np.arcsin(vector[2] / np.linalg.norm(vector))
theta = np.arccos(vector[0] / np.linalg.norm(xyproj))
if vector[1] < 0:
theta *= -1
trans = np.dot(self.rot_y(phi), self.rot_z(-theta))
return trans
[docs]
def align_plane(self, vector):
"""
Rotation around x axis to set a vector in the xy plane.
Parameters
==========
vector: array
vector to be rotated to be in the xy plane.
Return
======
(numpy.array) [3float x 3float] Transformation matrix.
"""
reference = vector.copy()
reference[0] = 0
angle = np.arccos(reference[1] / np.linalg.norm(reference))
if reference[2] < 0:
angle *= -1
return self.rot_x(-angle)
[docs]
def apply_trans(self, trans, indexes=None, shift=None):
"""
Apply a transformation to all vector positions of some atoms.
Parameters
==========
trans: array (3x3)
transformation matrix to be applied to all atom positions.
indexes: list or array. Defatult=None
indexes of the atoms to apply the transformation. Default None
that means the transformation is applied to the positions of all
the atoms.
shift: list. Default=None
translation vector added after the transformation. It is applied
to all the atoms, also to those not included in indexes.
Return
======
(list)[#natoms x 3float] new xyz positions of the N atoms. It
changes the positions of the internal atoms object.
"""
if shift is None:
shift = [0, 0, 0]
if indexes is None:
indexes = list(range(len(self.atoms)))
new_positions = []
for i, atom in enumerate(self.atoms):
if i in indexes:
new_positions.append(np.dot(trans, atom.position) + shift)
else:
new_positions.append(atom.position + shift)
self.atoms.set_positions(new_positions)
return new_positions
[docs]
def xy_alignment(self, index1, index2, index3=None, center=None):
"""
Transform the positions of the atoms such that the atoms of indexes1
and index2 are aligned in the x axis. The atom with index3 would be in
the xy plane in case to be given.
Parameters
==========
index1 and index2: int
indexes of the atoms to be aligned with the x-axis. The positive
direction of x would go from atom 1 to atom 2.
index3: int. Default=None
The atom with index 3 would be in the xy plane in case to be given.
center: int. Default=None
It must be index1 or index2, that means the atom with this index
will be placed in the origin. In case center=None (default), the
origin would be in the geometrical center between atoms with index1
and index2.
Return
======
(numpy.array)[#natoms x 3float] new xyz positions of the N atoms. It
changes the positions of the internal atoms object.
"""
if center == index1 or center == index2:
center = self.atoms[center].position
else:
pos1 = self.atoms[index1].position
pos2 = self.atoms[index2].position
center = (pos1 + pos2) / 2
self.atoms.set_positions(self.atoms.positions - center)
axis = self.atoms[index2].position - self.atoms[index1].position
self.apply_trans(self.align_axis(axis))
if index3 is not None:
third = self.atoms[index3].position
self.apply_trans(self.align_plane(third))
return self.atoms.positions
[docs]
def increase_distance(self, constraints, deltad):
"""
increase the distance between two atoms by moving them and keeping the
rest of the atoms in the same place.
Parameters
==========
constraints:
constraints with the shape (n, 2), where n is the number of
constraints and the first pair is the one to increase the
distance.
deltad: float
amount to add to the distance between atoms.
Return
======
(ase.Atoms) Internal Atoms object with the corresponding modification.
"""
d1norm = self.atoms.get_distance(constraints[0][0], constraints[0][1])
self.atoms.set_distance(constraints[0][0], constraints[0][1],
d1norm + deltad)
return self.atoms
[docs]
def increase_distance_with_constraints(self, constraints, deltad):
"""
Take a configuration and increase the distance between two atoms by
moving those atoms and the atoms connected to them through the other
constraints, and keeping the rest of the atoms in the same place. The
molecule is first centered and rotated such that the two atoms lie on
the x axis (see MoleculeSetter.xy_alignment).
Parameters
==========
constraints:
constraints with the shape (n, 2), where n is the number of
constraints and the first pair is the one to increase the
distance.
deltad: float
amount to add to the distance between atoms.
Return
======
(ase.Atoms) Internal Atoms object with the corresponding modification.
"""
self.xy_alignment(constraints[0][0], constraints[0][1])
left = [constraints[0][0]]
right = [constraints[0][1]]
len_right = 0
len_left = 0
while ((len_left != len(left)) and (len_right != len(right))):
len_left = len(left)
len_right = len(right)
for cons in constraints[1:]:
if ((cons[0] in left) and (cons[1] not in left)):
left.append(cons[1])
elif ((cons[1] in left) and (cons[0] not in left)):
left.append(cons[0])
elif ((cons[0] in right) and (cons[1] not in right)):
right.append(cons[1])
elif ((cons[1] in right) and (cons[0] not in right)):
right.append(cons[0])
print("These will be the atoms moved to the left: ", left)
print("These will be the atoms moved to the right: ", right)
new_positions = []
for i, atom in enumerate(self.atoms):
if i in right:
deltapos = np.array([deltad / 2, 0, 0])
new_positions.append(atom.position + deltapos)
elif i in left:
deltapos = np.array([-deltad / 2, 0, 0])
new_positions.append(atom.position + deltapos)
else:
new_positions.append(atom.position)
self.atoms.set_positions(new_positions)
return self.atoms
[docs]
def scale_distance(self, constraints, deltad, index3=None):
"""
Increase the distance between two atoms by aligning those atoms with
the x axis (the first atom at the origin) and scaling the x-coordinate
of all the atoms.
Parameters
==========
constraints: list
constraints with the shape (n, 2), where n is the number of
constraints and the first pair is the one to increase the
distance.
deltad: float
amount to add to the distance between atoms.
index3: int. Default=None
The atom with index 3 would be in the xy plane in case to be given.
Return
======
(list)[#natoms x 3float] scaled xyz positions, before the rotation of
index3 into the xy plane. The internal atoms object is also modified
(including that rotation).
"""
index1, index2 = constraints[0]
d1norm = self.atoms.get_distance(index1, index2)
self.xy_alignment(index1, index2, center=index1)
scale_factor = np.array([(d1norm + deltad) / d1norm, 1, 1])
new_positions = [atom.position * scale_factor
for atom in self.atoms]
self.atoms.set_positions(new_positions)
if index3 is not None:
third = self.atoms[index3].position
self.apply_trans(self.align_plane(third))
return new_positions
[docs]
class Alignment:
"""
Set of tools to apply the necessary transformations for aligning a the
three first eigenvectors of PCA with the x, y, z axis, in that order.
Parameters
==========
atoms: ase.Atoms
molecule to be aligned.
"""
def __init__(self, atoms):
"""
Stores the molecule to be aligned and a MoleculeSetter to transform it.
See the class docstring for the parameters.
"""
self.atoms = atoms
self.ms = MoleculeSetter(self.atoms)
[docs]
@staticmethod
def pca_vectors(atoms, indexes='all'):
"""
Find the eigen vectors in the Principal Component Analysis.
Parameters
==========
atoms: ase.Atoms
molecule to be analyzed.
indexes: Default="all"
subset of atoms to be included in the pca analysis. numbering
convention starting from 0.
Return
======
(array) principal components.
"""
try:
from sklearn.decomposition import PCA
except ModuleNotFoundError:
raise ModuleNotFoundError("This code requires sklearn")
if indexes == 'all':
indexes = list(range(len(atoms)))
positions = atoms.positions[indexes]
pca = PCA(n_components=3)
pca.fit(positions)
components = pca.components_
return components
[docs]
@staticmethod
def center_geo(atoms):
"""
Moves the coordinates origin to the geometrical center atoms.
Parameters
==========
atoms: ase.Atoms
Molecule to be moved.
Return
======
(ase.Atoms) Transformed atoms
"""
ms = MoleculeSetter(atoms)
geo_center = np.sum(atoms.positions, axis=0) / len(atoms)
ms.apply_trans(np.identity(3), shift=-geo_center)
return ms.atoms
[docs]
@staticmethod
def align_with_components(atoms):
"""
Aligns the coordinates x, y, z axis with the main vectors (in that
order).
Parameters
==========
atoms: ase.Atoms
Molecule to be aligned.
Return
======
(ase.Atoms) rotated and shifted molecule.
"""
ms = MoleculeSetter(atoms)
components = Alignment.pca_vectors(atoms)
ms.atoms += Atom('H', position=[0, 0, 0])
for vec_comp in components:
ms.atoms += Atom('H', position=vec_comp)
ms.xy_alignment(-4, -3, -2)
ms.atoms = ms.atoms[:-4]
return ms.atoms
[docs]
class PCAMatcher:
"""
Tools to align two molecules with their PCA main vectors.
Parameters
==========
reference: ase.Atoms
reference of the molecule to be matched.
to_compare: ase.Atoms
molecule to be compared to the reference.
"""
def __init__(self, reference, to_compare):
"""
Stores copies of the two molecules to be matched. They must have the same
number of atoms. See the class docstring for the parameters.
"""
assert len(reference) == len(to_compare), "The reference and the " + \
"structure to compare have to have the same number of atoms"
self.reference = reference.copy()
self.to_compare = to_compare.copy()
[docs]
def align_molecules(self):
"""
Align each molecule with their own main PCA vectors.
Return
======
(tuple) aligned reference and molecule to compare.
"""
self.reference = Alignment.center_geo(self.reference)
self.reference = Alignment.align_with_components(self.reference)
self.to_compare = Alignment.center_geo(self.to_compare)
self.to_compare = Alignment.align_with_components(self.to_compare)
return self.reference, self.to_compare
[docs]
def matching(self):
"""
Finds the order of the indices of the molecule which fit with the atoms
of the reference. So, if the outcome of this function is "out", then
PCAMatcher.to_compare[out[i]] == PCAMatcher.reference[i] should be
true.
Return
======
(numpy.array) indices of the molecule that match the indices of the
reference.
Note
----
Both molecules are aligned (see PCAMatcher.align_molecules) before
matching.
"""
self.align_molecules()
n_atoms = len(self.reference)
correspondence = -np.ones(n_atoms, dtype=int)
for i in range(n_atoms):
distances = np.linalg.norm(self.reference.positions[i]
- self.to_compare.positions,
axis=1)
i_min = np.where(distances == min(distances))[0][0]
if correspondence[i] == -1:
correspondence[i] = int(i_min)
if -1 in correspondence:
print('Warning: check repetitions because the next atoms of the '
'reference did not get a correspondent in the test:',
np.where(correspondence == -1)[0])
return correspondence
[docs]
def vmd_connectivity(mol):
"""
Extracts the connectivity of a molecule using VMD.
Parameters
==========
mol: str
path to a molecule file that VMD can read (e.g. xyz or pdb).
Return
======
(list) [n_atoms] list where the element i contains the indices of the
atoms bonded to the atom i. Indices start from 0 and every bond appears
twice (j in the list of i and i in the list of j).
Note
====
VMD must be available in the PATH. It uses 'eproducer extract_bonds_vmd'.
"""
raw = output_terminal(f'eproducer extract_bonds_vmd {mol} | ' +
r"awk '/atomselect0/{flag=1} /Info\)/{flag=0} flag' | tail -n +2",
print_output=False).split('\n')
connectivity = [[int(i) for i in line.split()] for line in raw[:-1]]
return connectivity
[docs]
def vmd_connectivity_unique(mol):
"""
Extracts the connectivity of a molecule using VMD, keeping each bond only
once.
Parameters
==========
mol: str
path to a molecule file that VMD can read (e.g. xyz or pdb).
Return
======
(list) [n_atoms] list where the element i contains the indices j > i of
the atoms bonded to the atom i. Indices start from 0.
"""
redundant = vmd_connectivity(mol)
unique = []
for i, connections in enumerate(redundant):
i_row = []
for j in connections:
if j > i:
i_row.append(j)
unique.append(i_row)
return unique
[docs]
def vmd_connectivity_matrix(mol):
"""
Builds the adjacency matrix of a molecule using the connectivity from VMD.
Parameters
==========
mol: str
path to a molecule file that VMD can read (e.g. xyz or pdb).
Return
======
(numpy.array) [n_atoms x n_atoms] symmetric matrix with 1 where two atoms
are bonded and 0 otherwise.
"""
redundant = vmd_connectivity(mol)
matrix = np.zeros((len(redundant), len(redundant)))
for i, connections in enumerate(redundant):
for j in connections:
matrix[i][j] = 1
return matrix