"""
Utilities for geometry manipulation or property calculation
Copyright Schrodinger, LLC. All rights reserved.
"""
import numpy
from scipy import optimize
from schrodinger import geometry as infrageom
from schrodinger.application.matsci import graph
from schrodinger.application.matsci import msprops
from schrodinger.structutils import analyze
[docs]def get_backbone_atoms(struct, exclude_hydrogen=True, max_fix=True):
"""
Return dictionary of backbone atoms for each molecule in the structure.
:param `schrodinger.structure.Structure` struct: The structure to find
backbone atoms in.
:param bool exclude_hydrogen: If true hydrogen will be excluded in
creation of graph and searching for backbone atoms. This results in
much faster calculation.
:param max_fix: Whether to fix internal degeneracy when maximum number of
path checks are reached
:type max_fix: bool
:return dict: the key of the dictionary is molecule number and the
value is the list of backbone atoms. If there are less than two atoms
in the molecule the associated value will be an empty list.
"""
heavy_atom_asl = 'heavy_atoms'
backbone = {}
# Create graph
atom_list = None
if exclude_hydrogen:
atom_list = analyze.evaluate_asl(struct, heavy_atom_asl)
struct_graph = analyze.create_nx_graph(struct, atoms=atom_list)
# Iterate through molecule to find backbone in each
for sub_graph in graph.get_sub_graphs(struct_graph):
# Get molecule number
nodes = list(sub_graph.nodes())
rep_node = nodes[0]
mol_num = struct.atom[rep_node].molecule_number
# Find backbone path
backbone_path = graph.find_backbone(sub_graph, max_fix=max_fix)
# Add to the dictionary if path is found, for single atoms there will be
# no path
if backbone_path:
backbone[mol_num] = backbone_path
return backbone
[docs]def get_ordered_polymer_backbone(struct, backbone_path, remove_side_chain):
"""
Backbone of a polymer chain should follow H--TH--TH--T format, where
backbone starts with head and ends with tail.
:param `schrodinger.structure.Structure` struct: The structure to find
backbone atoms in.
:param list backbone_path: list of atom indexes in the backbone
:param bool remove_side_chain: If true it will remove the side chain atoms
from the backbone in the start and end.
:return list: list of atom indexes in the backbone such that first atom is
always a head.
"""
ht_prop_name = [msprops.HEAD_ATOM_PROP, msprops.TAIL_ATOM_PROP]
# Check if head and tail information is present in backbone atoms
for atom_id in backbone_path:
# Check if atom is head or tail
atom = struct.atom[atom_id]
is_head, is_tail = [atom.property.get(x) for x in ht_prop_name]
if is_head or is_tail:
break
else:
# return unchanged list if no head and tail are available
return backbone_path
# In case tail is first backbone needs to be reversed
if is_tail:
backbone_path.reverse()
if remove_side_chain:
# Find index of first head atom and last tail atom in the backbone
ht_index = {}
found_prop, not_found_prop = [
ht_prop_name[0] if x else ht_prop_name[1]
for x in [is_head, is_tail]
]
ht_index[found_prop] = backbone_path.index(atom_id)
for index, atom_id in enumerate(backbone_path):
if struct.atom[atom_id].property.get(not_found_prop):
ht_index[not_found_prop] = index
if is_tail:
# Break only if finding first head
break
# Remove initial and final side chain atoms by selecting atoms between
# first head atom and last tail atom in the backbones
backbone_path = backbone_path[
ht_index[ht_prop_name[0]]:ht_index[ht_prop_name[1]]]
return backbone_path
[docs]def get_ordered_backbone_atoms(struct):
"""
Return dictionary of backbone atoms for each molecule in the structure.
:param `schrodinger.structure.Structure` struct: The structure to find
backbone atoms in.
:return dict(list): the key of the dictionary is molecule number and the
value is the list of backbone atoms
"""
fixed_struct = struct.copy()
# Get backbone atoms
backbone = get_backbone_atoms(fixed_struct)
# Fix order for polymeric systems
for mol_num, backbone_path in backbone.items():
backbone[mol_num] = get_ordered_polymer_backbone(
fixed_struct, backbone_path, False)
return backbone
[docs]def radius_of_gyration(struct, molecule_coms=False):
"""
Calculate radius of gyration for a structure
:param `structure.Structure` struct: The structure to calculate RG for
:param bool molecule_coms: Whether to use center of mass of each molecule in
the structure for calculating Rg
:rtype: float
:return: The radius of gyration
"""
xyzs = struct.getXYZ()
weights = [atom.atomic_weight for atom in struct.atom]
center_of_mass = numpy.average(xyzs, weights=weights, axis=0)
if molecule_coms:
# Method from https://pubs.acs.org/doi/abs/10.1021/acs.jpcb.0c02856
mol_coms = []
for mol in struct.molecule:
mol_xyzs = [atom.xyz for atom in mol.atom]
atomic_weights = [atom.atomic_weight for atom in mol.atom]
mol_coms.append(
numpy.average(mol_xyzs, weights=atomic_weights, axis=0))
mol_coms = numpy.array(mol_coms)
sum_distance_squared = numpy.sum(numpy.square(mol_coms -
center_of_mass))
r_g = numpy.sqrt(sum_distance_squared / len(struct.molecule))
else:
distance_squared = numpy.sum(numpy.square(xyzs - center_of_mass),
axis=1)
r_g = numpy.sqrt(
numpy.sum(weights * distance_squared) / struct.total_weight)
return r_g
[docs]def fit_structure_to_ellipsoid(struct):
"""
Fit the structure to an ellipsoid and return the ellipsoid semi-axes and
principal moments of inertia
The method is taken from https://pubs.acs.org/doi/abs/10.1021/jp047138e
and uses the following equations to calculate the semi-axes values:
I1 = M/5 * (a**2 + b**2)
I2 = M/5 * (a**2 + c**2)
I3 = M/5 * (b**2 + c**2)
Also see https://en.wikipedia.org/wiki/Ellipsoid#Dynamical_properties
:param `structure.Structure` struct: The structure to fit to an ellipsoid
:rtype: tuple of (float, float, float, [float, float, float])
:return: Tuple of semi-axes values in Angstroms and principal moments of
inertia in u * A**2
"""
moments = infrageom.principal_moments(struct).moments
moments = sorted(moments, reverse=True)
multiplier = struct.total_weight / 5
coefficients = numpy.array([[1, 1, 0], [1, 0, 1], [0, 1, 1]]) * multiplier
vals_squared = numpy.linalg.solve(coefficients, moments)
# If the atoms are in a plane, c^2 might be -1e-16 or something
vals_squared[vals_squared < 0] = 0
a_val, b_val, c_val = numpy.sqrt(vals_squared)
return a_val, b_val, c_val, moments
[docs]def fit_points_to_circle(x_vals, y_vals):
"""
Fit 2D points to a circle
:param numpy.array x_vals: Array containing point x values
:param numpy.array y_vals: Array containing point y values
:rtype: tuple(float, float, float)
:return: Center x value, center y value, and radius
"""
def calc_radii(x_c, y_c):
"""
calculate the distance of each point from the center (x_c, y_c)
"""
return numpy.sqrt((x_vals - x_c)**2 + (y_vals - y_c)**2)
def func_to_minimize(args):
"""
calculate the distance between the points and the mean circle centered
at (x_c, y_c)
"""
radii = calc_radii(*args)
return radii - radii.mean()
initial_guess = (0, 0)
results, _ = optimize.leastsq(func_to_minimize, initial_guess)
x_c, y_c = results
radius = calc_radii(*results).mean()
return x_c, y_c, radius