"""
Contains classes for working with Z-Matrices
Copyright Schrodinger, LLC. All rights reserved.
"""
# Contributors: Dave Giesen
from copy import deepcopy
import numpy
from schrodinger import structure
from schrodinger.application.jaguar import input as jaguar_input
from schrodinger.application.matsci import pbcutils
from schrodinger.infra import mm
from schrodinger.structutils.transform import get_normalized_vector
[docs]def rotation_matrix(axis, angle):
"""
Get rotation matrix based on the Euler-Rodrigues formula.
:param list axis: Axis, defined by 3 floats, will be normalized
:param float angle: Angle (deg) to rotate around the axis
:rtype: numpy.array
:return: Rotation matrix (3 x 3)
"""
# See: https://en.wikipedia.org/wiki/Euler%E2%80%93Rodrigues_formula
# This matrix is different compared to one computed with
# structutils.transform.get_rotation_matrix
# It seems that the difference comes from the left/right handed rotation
# Normalize the axis
axis = get_normalized_vector(numpy.array(axis))
angle = numpy.deg2rad(angle)
a = numpy.cos(angle / 2)
b, c, d = -axis * numpy.sin(angle / 2)
a2, b2, c2, d2 = numpy.power([a, b, c, d], 2)
return numpy.array(
[[a2 + b2 - c2 - d2, 2 * (b * c - a * d), 2 * (b * d + a * c)],
[2 * (b * c + a * d), a2 + c2 - b2 - d2, 2 * (c * d - a * b)],
[2 * (b * d - a * c), 2 * (c * d + a * b), a2 + d2 - b2 - c2]]
) # yapf: disable
[docs]class ZAtom(object):
"""
A single atom in a Z-matrix
"""
COORD_LETTERS = ['r', 'a', 'd']
[docs] def __init__(self, index, element, values, constraints, neighbor_indexes,
neighbor_elements):
"""
Initialize ZAtom object
:param int index: Atom index with the structure
:param str element: Atom element
:param list values: Atom values (bond distance, angle, torsion)
:param set constraints: Atom constraints
:param list neighbor_indexes: Atom neighboring indices
:param list neighbor_elements: Atom neighboring elements
"""
self.index = index
self.element = element
self.values = list(values)
self.constraints = set(constraints)
self.neighbor_indexes = list(neighbor_indexes)
self.neighbor_elements = list(neighbor_elements)
[docs] @classmethod
def fromJagin(cls, jagin, index):
"""
Create a ZAtom object from jaguar input.
:type jagin: `schrodinger.application.jaguar.input.JaguarInput`
:param jagin: The JaguarInput object the Z-matrix is for
:type index: int
:param index: The atom index (1-based) that this ZAtom is for
:rtype: ZAtom
:return: New ZAtom object
"""
struct = jagin.getStructure()
atom = struct.atom[index]
ZMAT_NUM = 0
(ctype, ind1, val1, ind2, val2, ind3,
val3) = mm.mmjag_zmat_atom_get(jagin.handle, ZMAT_NUM, atom.index)
neighbor_indexes = [x for x in [ind1, ind2, ind3] if x]
neighbor_elements = [struct.atom[x].element for x in neighbor_indexes]
# The value for an unused coordinate is 0, but 0 is also a possible
# value for real angles and dihedrals. Zipping the values and indexes
# together truncates the value list at the number of real coordinates.
valdexes = list(zip([val1, val2, val3], neighbor_indexes))
values = [x for x, y in valdexes]
index = atom.index
element = atom.element
constraints = set()
return cls(index, element, values, constraints, neighbor_indexes,
neighbor_elements)
def __deepcopy__(self, memo):
"""
Deep copy a ZAtom.
:param memo: Memo to be passed to the deepcopy calls (if any)
:type memo: Dict or None
:rtype: ZAtom
:return: Newly copied ZAtom
"""
cls = type(self)
newobj = cls(self.index, self.element, self.values, self.constraints,
self.neighbor_indexes, self.neighbor_elements)
return newobj
[docs] def constrainCoord(self, coord):
"""
Constrain the given coord so it can't be optimized
:type coord: int
:param coord: 0 (bond), 1 (angle) or 2 (dihedral)
"""
self.constraints.add(coord)
[docs] def freeCoord(self, coord):
"""
Allow the given coord to be optimized
:type coord: int
:param coord: 0 (bond), 1 (angle) or 2 (dihedral)
"""
self.constraints.remove(coord)
[docs] def constrainAll(self):
"""
Constrain all the coords so they can't be optimized
"""
self.constraints.update(range(len(self.COORD_LETTERS)))
[docs] def freeAll(self):
"""
Allow all the coords to be optimized
"""
self.constraints.clear()
[docs] def bumpIndexes(self, bump):
"""
Modify all the indexes (both for self and any referenced atoms) by a
constant number. Used when this z-matrix is being added to another, and
now atom 1 is actually atom 20...
:type bump: int
:param bump: The amount to modify each index
"""
self.index = self.index + bump
self.neighbor_indexes = [x + bump for x in self.neighbor_indexes]
[docs] def getSymbolicLine(self):
"""
Get the line that defines the atoms for the bond, angle and dihedral and
also the coordinate names. i.e.
C11 H1 r11 C2 a11 C7 d11
:rtype: str
:return: The string defining the z-matrix coordinates for this atom.
Includes the trailing return character
"""
tokens = [self.getAtomName(self.element, self.index)]
for coord, (element, index) in enumerate(
zip(self.neighbor_elements, self.neighbor_indexes)):
neighbor_name = self.getAtomName(element, index)
coord_name = self.getCoordName(coord)
tokens.extend([neighbor_name, coord_name])
return ' '.join(tokens) + '\n'
[docs] def getVariableLines(self):
"""
Get a series of lines that defines each coordinate, i.e.::
r11 = 1.452
a11 = 102.4
d11 = 180
:rtype: str
:return: The string defining the z-matrix coordinate values for this
atom. Includes the trailing return character at the end of the block.
"""
lines = ""
for index, value in enumerate(self.values):
cname = self.getCoordName(index)
if index in self.constraints:
hashtag = '#'
else:
hashtag = ""
lines += cname + ' = ' + "%.12f" % value + hashtag + '\n'
return lines
[docs] def getAtomName(self, element, index):
"""
Get the name of an atom with the given element and index
:type element: str
:param element: The atomic symbol of the atom
:type index: int
:param index: The index of the atom
:rtype: str
:return: The name of the atom, such as "C12"
"""
return element + str(index)
[docs] def getCoordName(self, coord):
"""
Get the name of a coordinate for this atom, such as r11 or d7
:type coord: int
:param coord: 0 (bond), 1 (angle) or 2 (dihedral)
:rtype: str
:return: The name of the coordinate, such as "r12"
"""
return self.COORD_LETTERS[coord] + str(self.index)
[docs] def addCoord(self, index, element, value):
"""
Add an internal coordinate for this atom referenced to a target atom.
Note that coordinates must be added in the order bond, angle, dihedral.
:type index: int
:param index: The atom index of the target atom (such as the atom the
bond extends to)
:type element: str
:param element: The element of the target atom
:type value: float
:param value: The value of this coordinate
"""
self.neighbor_indexes.append(index)
self.neighbor_elements.append(element)
self.values.append(value)
[docs] def changeCoord(self, coord, value, index=None, element=None):
"""
Change the definition of an existing internal coordinate. If index and
element are passed in the target atom for the coordinate will be
changed.
:type coord: int
:param coord: 0 (bond), 1 (angle) or 2 (dihedral)
:type value: float
:param value: The value of this coordinate
:type index: int
:param index: The atom index of the target atom (such as the atom the
bond extends to), 1-based. Both index and element must be supplied for
them to be used.
:type element: str
:param element: The element of the target atom. Both index and element
must be supplied for them to be used.
"""
if index and element:
self.neighbor_indexes[coord] = index
self.neighbor_elements[coord] = element
self.values[coord] = value
[docs]class ZMatrix(object):
"""
Contains the z-matrix for a structure
"""
[docs] def __init__(self, struct=None):
"""
Create a ZMatrix instance
:type struct: `schrodinger.structure.Structure` or None
:param struct: The structure this zmatrix is for
"""
if struct:
jagin = jaguar_input.JaguarInput(structure=struct, name='temp')
jagin.makeInternalCoords()
self.atoms = [
ZAtom.fromJagin(jagin, x)
for x in range(1,
jagin.getAtomCount() + 1)
]
self.setBasisVectors(struct)
else:
self.atoms = []
[docs] def setBasisVectors(self, struct):
"""
Set origin and basis vectors. These can be used to recover original
Cartesian coordinates from Z-matrix.
:param structure.Structure: Input structure
"""
# See getStructure function on how these are used
if struct.atom_total >= 3:
self.origin = numpy.array(struct.atom[1].xyz)
self.bond1_vec = get_normalized_vector(
numpy.asarray(struct.atom[2].xyz) - self.origin)
self.bond2_vec = get_normalized_vector(
numpy.asarray(struct.atom[3].xyz) -
numpy.asarray(struct.atom[2].xyz))
else:
self.origin = numpy.zeros(3) # Origin at zero
self.bond1_vec = numpy.array([1, 0, 0]) # First vector along X
self.bond2_vec = numpy.array([0, 1, 0]) # Second vector along Y
[docs] def getStructure(self):
"""
Create and return structure based on the zmatrix.
:rtype: `schrodinger.structure.Structure`
:return: Newly created structure
"""
struct = structure.create_new_structure()
for idx, zatom in enumerate(self.atoms):
# Note that indices can appear not sorted ([1, 0] vs [0, 1])
indices = numpy.array(zatom.neighbor_indexes) - 1
if idx == 0:
# First atom, place at the origin
struct.addAtom(zatom.element, *self.origin)
continue
elif idx == 1:
# Second atom, place along the bond1_vec vector
dist, = zatom.values
vec = struct.getXYZ()[0] + dist * self.bond1_vec
struct.addAtom(zatom.element, *vec)
continue
elif idx == 2:
# Third atom, place along the bond vector and angle
dist, angle = zatom.values
a1_vec, a2_vec = struct.getXYZ()[indices]
# Vector has origin at neighbor atom and direction along the
# bond2_vec
vec = a1_vec + dist * self.bond2_vec
atom = struct.addAtom(zatom.element, *vec)
# Adjust angle, atom indices that define it are reversed
atom_indices = list(reversed(indices + 1))
struct.adjust(angle, atom_indices[0], atom_indices[1],
atom.index)
continue
# Forth atom and beyond
dist, angle, dihedral = zatom.values
vec1, vec2, vec3 = struct.getXYZ()[indices]
vec_a, vec_b = vec2 - vec1, vec2 - vec3
dist *= get_normalized_vector(vec_a)
normal = numpy.cross(vec_a, vec_b)
tmatrix1 = rotation_matrix(normal, angle)
dist = numpy.dot(tmatrix1, dist)
tmatrix2 = rotation_matrix(vec_a, dihedral)
dist = numpy.dot(tmatrix2, dist)
vec = vec1 + dist
struct.addAtom(zatom.element, *vec)
return struct
[docs] def fromStructure(self, struct):
"""
Given the input structure, generate a ZMatrix with the same atom
neighbors lists as self.
:param structure.Structure struct: The structure this zmatrix is for
:rtype: ZMatrix
:return: Newly created ZMatrix
"""
# Note that sometimes for two structures with the same list of atoms
# but different coordinates, different zmatrices are generated by Jaguar
# hence this function
assert ([a.element for a in self.atoms
] == [a.element for a in struct.atom])
# Z-Matrix built using makeInternalCoords() in the constructor is not
# PBC aware. To match, remove PBC properties here
tmp_st = struct.copy()
pbcutils.remove_pbc_props(tmp_st)
newobj = type(self)()
for atom in self.atoms:
new_atom = deepcopy(atom)
for idx, nidx in enumerate(atom.neighbor_indexes):
new_atom.values[idx] = tmp_st.measure(
atom.index, *atom.neighbor_indexes[:idx + 1])
newobj.atoms.append(new_atom)
newobj.setBasisVectors(struct)
return newobj
[docs] def copy(self):
"""
Copy Z-matrix.
:rtype: ZMAtrix
:return: Copied Z-matrix
"""
newobj = type(self)()
newobj.atoms = [deepcopy(atom) for atom in self.atoms]
newobj.origin = numpy.array(self.origin)
newobj.bond1_vec = numpy.array(self.bond1_vec)
newobj.bond2_vec = numpy.array(self.bond2_vec)
return newobj
[docs] def difference(self, other_zmat):
"""
Compute difference Z-matrix between other zmatrix and self zmatrix
(in this order).
:param ZMatrix other_zmat: Other Z-Matrix
:rtype: ZMatrix
:return: Difference Z-matrix
"""
def diff_angles(angles1, angles2):
"""
Get the smallest angle difference between -180 and 180.
Pairs of angles are defined as angle2 - angle1.
"""
# Example. For y, x = 359, 2
# result is -3
# See also https://stackoverflow.com/a/7869457
return [(y - x + 180) % 360 - 180 for x, y in zip(angles1, angles2)]
assert len(self.atoms) == len(other_zmat.atoms)
ret = ZMatrix()
for atom1, atom2 in zip(self.atoms, other_zmat.atoms):
new_atom = deepcopy(atom1)
if len(atom1.values):
# at least bond is present
new_atom.values[0] = atom2.values[0] - atom1.values[0]
if len(new_atom.values) >= 2:
new_atom.values[1:] = diff_angles(atom1.values[1:],
atom2.values[1:])
ret.atoms.append(new_atom)
ret.origin = other_zmat.origin - self.origin
ret.bond1_vec = other_zmat.bond1_vec - other_zmat.bond1_vec
ret.bond2_vec = other_zmat.bond2_vec - other_zmat.bond2_vec
return ret
[docs] def interpolate(self, other_zmat, disps, indices=None):
"""
Interpolate between two Z-matrices. Initial Z-matrix is self, other
zmatrix is final point.
:param ZMatrix other_zmat: Final point of zmatrix
:param list disps: List of displacement coefficients
:type indices: list[int] or None
:param indices: List of atom indices to move, if None, will move all
atoms
:rtype: list[ZMatrix]
:return: List of interpolated Z-matrices
"""
assert len(self.atoms) == len(other_zmat.atoms)
indices = set(
range(1,
len(self.atoms) + 1) if indices is None else indices)
indices = sorted(indices)
zmat_diff = self.difference(other_zmat)
zmats = []
for disp in disps:
zmat = self.copy()
for atom, atom_diff in zip(zmat.atoms, zmat_diff.atoms):
if atom.index in indices:
atom.values = (numpy.asarray(atom.values) +
numpy.asarray(atom_diff.values) * disp)
zmat.origin += zmat_diff.origin * disp
zmat.bond1_vec = get_normalized_vector(zmat.bond1_vec +
zmat_diff.bond1_vec * disp)
zmat.bond2_vec = get_normalized_vector(zmat.bond2_vec +
zmat_diff.bond2_vec * disp)
zmats.append(zmat)
return zmats
[docs] def bumpIndexes(self, bump):
"""
Bump all atom indexes for this z-matrix
:type bump: int
:param bump: The amount to modify each index
"""
for atom in self.atoms:
atom.bumpIndexes(bump)
[docs] def constrainAll(self):
"""
Constrain all internal coordinates
"""
for atom in self.atoms:
atom.constrainAll()
[docs] def freeAtom(self, index):
"""
Free all the internal coordinates for a single atom
:type index: int
:param index: The atom index (1-based) of the atom to free
"""
self.atoms[index - 1].freeAll()
[docs] def freeAtomBond(self, index):
"""
Free all the bond internal coordinate for a single atom
:type index: int
:param index: The atom index (1-based) of the atom to free
"""
self.atoms[index - 1].freeCoord(0)
[docs] def freeAtomAngle(self, index):
"""
Free all the angle internal coordinate for a single atom
:type index: int
:param index: The atom index (1-based) of the atom to free
"""
self.atoms[index - 1].freeCoord(1)
[docs] def freeAtomDihedral(self, index):
"""
Free all the angle internal coordinate for a single atom
:type index: int
:param index: The atom index (1-based) of the atom to free
"""
self.atoms[index - 1].freeCoord(2)
[docs] def findNonlinearAngle(self, struct, atom_a, atom_b, do_not_use,
index_less_than, linear):
"""
Find an angle between atom_a, atom_b and some other atom that is smaller
than linear
:type struct: `schrodinger.structure.Structure`
:param struct: The struct to use
:type atom_a: int
:param atom_a: The 1-indexed index of atom A for the A-B-C angle
:type atom_b: int
:param atom_b: The 1-indexed index of atom B for the A-B-C angle
:type do_not_use: set
:param do_not_use: An set of atom indexes that should not be used for
the C atom for the A-B-C angle
:type index_less_than: int
:param index_less_than: Only look for indexes less than this number
:type linear: float
:param linear: The threshold where any angle greater than this is
considered linear
:rtype: int, float or None, None
:return: The index of an atom that forms an angle with atom_a and atom_b
that is less than linear degrees, and the angle it forms. None, None
is returned if none of the indexes in choose_from gave a non-linear
angle.
"""
for index in range(1, index_less_than):
if index not in do_not_use:
angle = struct.measure(atom_a, atom_b, index)
if abs(angle) < linear:
return index, angle
return None, None
[docs] def merge(self, child_zmat, struct, glue_index):
"""
Combine the atom lists from this (parent) z-matrix and another (child)
z-matrix. This involves not just combining the lists, but also defining
the following new internal coordinates:
child_atom_1: bond, angle, dihedral
child_atom_2: angle, dihedral
child_atom_3: dihedral
All new internal coordinates will reference an atom from the parent
structure. These six new internal coordinates fully specify the
parent-child orientation.
:type child_zmat: `ZMatrix`
:param child_zmat: The child z-matrix to merge with this one
:type struct: `schrodinger.structure.Structure`
:param struct: The structure that results from merging these two
z-matrices. This is needed to compute the new bonds, angles and
dihedrals.
:type glue_index: int
:param glue_index: The index of the atom in the parent structure that
should be bonded to the first atom of the child structure. This index
should be 1-based (i.e. structure.atom indexed, not list indexed)
"""
parent_atom_total = len(self.atoms)
child_zmat.bumpIndexes(parent_atom_total)
# Compile a list of parent atoms to use to define the new child-parent
# orientation coordinates
parent_coord_indexes = [glue_index]
neighbors = self.atoms[glue_index - 1].neighbor_indexes[:]
# By preference use the atoms that the glue atom uses for coordinates
parent_coord_indexes.extend(neighbors)
if len(parent_coord_indexes) < 3:
# glue_index atom must be one of the first three, so it doesn't have
# three neighbors defined in the z-matrix - add the remaining atoms
# from the first three parent atoms (if they exist)
maxatoms = min(3, parent_atom_total)
for index in range(1, maxatoms + 1):
if index not in parent_coord_indexes:
parent_coord_indexes.append(index)
# Now add z-matrix coordinates to fill out the first 3 child atoms. Each
# new z-matrix coordinate will be referenced to an atom in the parent
# structure.
# LINEAR = Jaguar cutoff for rejecting linear angle coordinates
LINEAR = 165.
atoms_to_fix = min(3, len(child_zmat.atoms))
for child_index in range(atoms_to_fix):
# Note - child index is 0-indexed
child_atom = child_zmat.atoms[child_index]
# Note - child_atom.index is 1-indexed
coord_indexes = parent_coord_indexes[:]
if child_index < 1:
# Bond only for the first child atom
parent_index = coord_indexes.pop(0)
bond = struct.measure(child_atom.index, parent_index)
element = struct.atom[parent_index].element
child_atom.addCoord(parent_index, element, bond)
if child_index < 2 and coord_indexes:
# Angle only for the first and second child atoms
parent_index = coord_indexes.pop(0)
bond_index = child_atom.neighbor_indexes[0]
angle = struct.measure(child_atom.index, bond_index,
parent_index)
if abs(angle) > LINEAR:
# Angle A-B-C is too close to linear. We will definitely
# use A-B (they are bonded), so find a new C so that A-B-C
# doesn't give 180
do_not_use = set(parent_coord_indexes)
index, new_angle = self.findNonlinearAngle(
struct, child_atom.index, bond_index, do_not_use,
parent_atom_total + 1, LINEAR)
if index is not None:
c_index = index
abc_angle = new_angle
else:
c_index = parent_index
abc_angle = angle
a_index = child_atom.index
b_index = bond_index
print('Warning: Probe-Target orientation involves an '
'angle near 180, formed by %d, %d, %d atoms. '
'Geometry may be unstable.' %
(a_index, b_index, c_index))
else:
c_index = parent_index
abc_angle = angle
element = struct.atom[c_index].element
child_atom.addCoord(c_index, element, abc_angle)
if coord_indexes:
# Dihedral for all 3 first child atoms
# There are two angles involved in an A-B-C-D dihedral: A-B-C
# and B-C-D. We need to ensure that neither of those angles are
# greater than the Jaguar LINEAR cutoff. In the process of
# defining the angle, we already ensured that A-B-C was OK, so
# now we need to concern ourselves with B-C-D. We can ensure
# that B-C-D will be OK by using the atoms that B defines
d_index = coord_indexes.pop(0)
bcd_angle = struct.measure(child_atom.neighbor_indexes[0],
child_atom.neighbor_indexes[1],
d_index)
if abs(bcd_angle) > LINEAR:
do_not_use = set(parent_coord_indexes +
child_atom.neighbor_indexes)
index, angle = self.findNonlinearAngle(
struct, child_atom.neighbor_indexes[0],
child_atom.neighbor_indexes[1], do_not_use,
parent_atom_total + 1, LINEAR)
if index is None:
a_index = child_atom.neighbor_indexes[0]
b_index = child_atom.neighbor_indexes[1]
print('Warning: Probe-Target orientation involves a '
'torsion with one interior angle near 180, '
'formed by %d, %d, %d atoms. Geometry may be '
'unstable.' % (a_index, b_index, d_index))
else:
d_index = index
torsion = struct.measure(child_atom.index,
child_atom.neighbor_indexes[0],
child_atom.neighbor_indexes[1],
d_index)
element = struct.atom[d_index].element
child_atom.addCoord(d_index, element, torsion)
self.atoms = self.atoms + child_zmat.atoms
[docs]def add_zmatrix_lines(lines, zmat):
"""
Add the lines defining a z-matrix to the character string using the format
required by Jaguar input files.
:type lines: str
:param lines: The string to add the z-matrix to
:type zmat: `schrodinger.application.matsci.zmutils.ZMatrix`
:param zmat: The ZMatrix that defines the lines to be added
:rtype: str
:return: The original lines with the z-matrix lines appended. The initial
'&zmat' and the closing '&' are appended by this function
"""
lines += '&zmat\n'
for atom in zmat.atoms:
lines += atom.getSymbolicLine()
lines += '&\n&zvar\n'
for atom in zmat.atoms:
lines += atom.getVariableLines()
lines += '&'
return lines
[docs]def replace_coords_with_zmat(path, zmat):
"""
Replace the coordinate section in a jaguar input file with the z-matrix.
Replaces the current:
&zmat
...
&
section with
&zmat
...
&
&zvar
...
&
:type path: str
:param path: The path to the input file
@zmat: ZMatrix
:param: The ZMatrix object to use to replace the coordinates
"""
# Now read in the input file, replace the coordinate section with our
# z-matrix, and write it back out.
input_file = open(path, 'r')
all_lines = ""
skipping = False
for line in input_file:
if not line.startswith('&zmat'):
if not skipping:
# If we're currently outside the &zmat section add this line
all_lines += line
else:
if line.startswith('&'):
# We're inside the zmat section and just found the end of it
skipping = False
else:
# We just found the beginning of the zmat section
skipping = True
all_lines = add_zmatrix_lines(all_lines, zmat)
input_file.close()
with open(path, 'w') as input_file:
input_file.write(all_lines)