"""
Module to do geometry based structure calculation which works both on structure
and trajectory.
There are three levels. Base level is MSAtom level is single atom or group
atoms treated as as a single atom. This is similar to Coarse grain bead which
can represent multiple atoms.
Second level is MSAtomGroup which is a holds multiple atom group. This is the
class which is responsible for geometric calculation using MSAtom coordinates.
Third level is MSAtomGroups which holds multiple atom groups and runs the same
geometric calculation on all of them.
Copyright Schrodinger, LLC. All rights reserved.
"""
from dataclasses import dataclass, fields
from typing import Any
import numpy as np
from scipy import optimize
from scipy.stats import trim_mean
from schrodinger.application.matsci import graph
from schrodinger.application.matsci import msutils
from schrodinger.application.matsci import codeutils
from schrodinger.application.matsci.coarsegrain import FusedRing
from schrodinger.structutils import analyze
from schrodinger.structutils import transform
topo = codeutils.get_safe_package('desmond.topo')
analysis = codeutils.get_safe_package('desmond.analysis')
GROUP_NONE = 'atom'
GROUP_MASS = 'com'
GROUP_CENTROID = 'centroid'
HELIX_AIDS_DELIMITER = ';'
[docs]class MSAtom:
"""
A atom like object which holds an atom or group of atoms that can be
construed as a single atom.
"""
[docs] def __init__(self, indices, struct, pos_type=GROUP_NONE):
"""
Initialize MSAtom
:type index: int
:param index: index of the atom
:type struct: `schrodinger.structure.Structure` or `cms.Cms`
:param struct: Structure to which the atom belongs. It will be cms
model if trajectory is used, and structure if trajectory is not
used.
:type pos_type: str
:param pos_type: Select GROUP_MASS or GROUP_CENTROID group all the
atoms into single atom group represented by the center of mass and
geometric center respectively. GROUP_NONE will not group any of the
atom, however, indices should contain a single atom. Defaults to
GROUP_NONE.
:raises ValueError: If there are more than one atom with pos_type as
GROUP_NONE
"""
self.struct = struct
self.indices = indices
self.pos_type = pos_type
if len(self.indices) != 1 and self.pos_type == GROUP_NONE:
raise ValueError(
'MSAtom cannot be created with multiple atoms while using '
'no grouping.')
try:
self.gids = topo.aids2gids(self.struct, [self.indices],
include_pseudoatoms=False)[0]
except (AttributeError, ImportError):
self.gids = None
def _getAllXYZ(self, pos=None):
"""
Get coordinate values for all the indices in the atom
:type pos: `numpy.ndarray` or None
:param pos: all coordinates in the trajectory for the current/analyzed
frame. None if it is a static structure
:rtype: `numpy.ndarray`
:return: coordinate of all the atoms making this MSAtom
:raises AttributeError: Raise error when gids were not set and
coordinates were requested from the trajectory coordinates.
"""
if pos is not None:
if self.gids is None:
raise AttributeError('GID not set for the atom.')
return pos[self.gids]
else:
return np.array([self.struct.atom[x].xyz for x in self.indices])
[docs] def getXYZ(self, all_pos=None):
"""
Get current coordinate of the atom group
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
:rtype: `numpy.ndarray`
:return: current coordinate for the atom group
:raises AttributeError: Raise attribute error if grouping is None and
single xyz is requested
"""
if self.pos_type == GROUP_CENTROID:
return np.mean(self._getAllXYZ(all_pos), axis=0).tolist()
elif self.pos_type == GROUP_MASS:
masses = [self.struct.atom[x].atomic_weight for x in self.indices]
return np.average(self._getAllXYZ(all_pos), axis=0,
weights=masses).tolist()
else:
# Return the first set of coordinate from all coordinates when
# there is not grouping
return self._getAllXYZ(all_pos)[0]
@property
def index(self):
"""
Index of the atom in the structure if MSAtom has only one atom. Raises
error if there are multiple atoms, ie when using pos type
:rtype: int
:return: index of atom in the structure
:raises ValueError: If there are more than one atom making this MSAtom
"""
if len(self.indices) != 1:
raise ValueError(
'Index property not set for MSAtom containing more than one '
'atom.')
return self.indices[0]
@property
def atomic_weight(self):
"""
Get atomic weight of the atom
:rtype: float
:return: atomic weight of the atom
"""
return sum(self.struct.atom[x].atomic_weight for x in self.indices)
[docs]class HelixArc(MSAtom):
"""
A MSAtom that constitutes of 4 consecutive atoms representing an helix arc
"""
[docs] def getArcCenter(self, pos):
"""
Fit the arc points to a sphere and return the center of the sphere
:type pos: `numpy.ndarray`
:param pos: coordinates of the arc points
:rtype: `numpy.ndarray`
:return: coordinate of the center
"""
def _vec_dist_score(center):
"""
Calculate the algebraic distance between the arc points and the
mean sphere centered
:param `numpy.ndarray` center: Estimated center for the sphere to
calculate all the distances from
:rtype: `numpy.ndarray`
:return: distance between the arc points and the mean sphere
centered
"""
least_sq_distances = np.linalg.norm(np.subtract(pos, center),
axis=1)
return least_sq_distances - least_sq_distances.mean()
center_estimate = np.mean(pos, axis=0)
result = optimize.least_squares(_vec_dist_score, center_estimate)
center = result.x
return np.array(center)
[docs] def calculateHelixProp(self, all_pos=None):
"""
Calculate the local helix properties using 4 points in the arc of the
helix
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
"""
# Get atom coordinates and vectors
pos = np.array(self._getAllXYZ(all_pos))
if len(pos) != 4:
raise ValueError('Helix arc should be a group of four atoms.')
vecs = np.diff(pos, axis=0)
# Find angle bisector of local helix
# BIOPOLYMERS VOL.5. PP.673-679(1967) Eq 11 and 12
bisects = [vecs[x] - vecs[x + 1] for x in range(2)]
# Make sure bisectors are in same direction, pointing outwards
# This occurs when 2nd or 3rd lie perfectly in the arc to form concave
if np.dot(*bisects) < 0:
# Create vector between first and last atom
end_vec = pos[0] - pos[-1]
# Flip the bisector which pointing towards end vector
arc_check = [np.dot(x, end_vec) < 0 for x in bisects]
try:
flip_index = arc_check.index(False)
except ValueError:
# In case both are not pointing towards the vector
# It is zig-zag and not an arc
self.is_arc = False
return
bisects[flip_index] *= -1
self.is_arc = True
# Find arc center
self.sphere_center = self.getArcCenter(pos)
self.center = pos.mean(axis=0)
# Find local radius
avg_radius = np.mean([
transform.get_vector_magnitude(x - self.sphere_center) for x in pos
])
self.local_diameter = 2 * avg_radius
# Find twist angle
# BIOPOLYMERS VOL.5. PP.673-679(1967) Eq 14, 15, and 16
bisect_angle = transform.get_angle_between_vectors(*bisects)
local_costheta = np.cos(bisect_angle)
self.twist = np.rad2deg(bisect_angle)
# Find local axis vector
local_axis = np.cross(*bisects)
# Not an arc if bisectors are on the same plane
if transform.get_vector_magnitude(local_axis) == 0.0:
self.is_arc = False
return
else:
self.local_axis = local_axis / np.linalg.norm(local_axis)
# Find height
center_to_pos_vec = [pos[x] - self.sphere_center for x in [1, 2]]
costheta_vals = [
np.cos(transform.get_angle_between_vectors(x, self.local_axis))
for x in center_to_pos_vec
]
distances = [
transform.get_vector_magnitude(x) * y
for x, y in zip(center_to_pos_vec, costheta_vals)
]
arc_height = np.abs(np.diff(distances))[0]
self.rise = arc_height
# Per turn
self.base_per_turn = 360 / self.twist
# Find pitch
local_angle = np.arccos(local_costheta)
self.pitch = arc_height * 2 * np.pi / local_angle
###############################################################################
# Atom Group Types - Constructed using MSAtom
###############################################################################
[docs]class MSAtomGroup:
"""
Class to calculate properties of a group of MSAtoms
"""
BASE_ATOM_TYPE = MSAtom
[docs] def __init__(self, struct, group_atom_ids, pos_type=GROUP_NONE):
"""
Initialize MSAtomGroup
:type struct: `schrodinger.structure.Structure` or `cms.Cms`
:param struct: Structure to which the atom group belongs. It will be
cms model if trajectory is used, and structure if trajectory is not
used.
:type group_atom_ids: list
:param group_atom_ids: list of atom ids to add to the group
:type pos_type: str
:param pos_type: Select GROUP_MASS or GROUP_CENTROID group all the
atoms into single-atom group represented by the center of mass and
geometric center, respectively. GROUP_NONE will not group any of
the atoms. However, indices should contain a single atom when using
GROUP_NONE. Defaults to GROUP_NONE.
"""
self.struct = struct
self.group_atom_ids = group_atom_ids
self.pos_type = pos_type
self.group_atoms = []
self.groupAtoms()
[docs] def groupAtoms(self):
"""
Create atom groups from group atom ids
"""
for atom_index in self.group_atom_ids:
self.group_atoms.append(
self.BASE_ATOM_TYPE([atom_index], self.struct, self.pos_type))
[docs] def getGroupXYZ(self, all_pos=None):
"""
Get xyz values for the atom group
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
:rtype: array
:return: xyz values for the atom group
"""
return [x.getXYZ(all_pos) for x in self.group_atoms]
[docs] def getCentroid(self, all_pos=None):
"""
Get centroid for the atom group
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
:rtype: `numpy.ndarray`
:return: centroid for the atom group, it is same as the single atom xyz
when there is no grouping.
"""
return np.mean(self.getGroupXYZ(all_pos), axis=0).tolist()
[docs] def getCenterOfMass(self, all_pos=None):
"""
Get center of mass for the atom group
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
:rtype: `numpy.ndarray`
:return: center of mass for the atom group, it is same as the single
atom xyz when there is no grouping.
"""
masses = [x.atomic_weight for x in self.group_atoms]
return np.average(self.getGroupXYZ(all_pos), axis=0,
weights=masses).tolist()
[docs]class MSRingAtomGroup(MSAtomGroup):
"""
Class containing group of MSAtoms where each atoms that are part of the
same ring are grouped into a single MSAtom
"""
[docs] def groupAtoms(self):
"""
Overwrite parent group to create group where MSAtom that are part of
the same ring are grouped into a single MSAtom
"""
group_atom_idx = set(self.group_atom_ids)
fused_rings = FusedRing.getFusedRings(self.struct)
new_group = []
# Add fused ring in the atom group as a single group
for fring in fused_rings:
# Check if atom group contains any atoms from the current ring
ring_idx = set(fring.getAtomIndices())
interset = group_atom_idx.intersection(ring_idx)
if interset:
# Remove the atom indices from the atom index tracking that
# are part of the fused ring
group_atom_idx = group_atom_idx.difference(interset)
# Add the fused ring to the new group
new_group.append(
self.BASE_ATOM_TYPE(list(ring_idx),
self.struct,
pos_type=GROUP_CENTROID))
# Add individual atoms that are not part of fused rings to the atom
# group
for atom_idx in group_atom_idx:
new_group.append(self.BASE_ATOM_TYPE([atom_idx], self.struct))
self.group_atoms = new_group.copy()
[docs]@dataclass
class HelixData:
"""
Data class to store helix data values
:type time: float
:param time: time associated with helix
:type diameter: float
:param diameter: diameter of the helix
:type twist: float
:param twist: twist angle of the helix
:type rise: float
:param rise: rise distance of each helix
:type base_per_turn: float
:param base_per_turn: average number of bases in each turn
:type pitch: float
:param pitch: pitch of the helix
:type narcs: int
:param narcs: number of arc used to define the helix
:type aids: str
:param aids: string of atom indices used to generate helix separated by
semi-colon
"""
time: Any = None
diameter: float = None
twist: float = None
rise: float = None
base_per_turn: float = None
pitch: float = None
narcs: int = None
aids: str = None
[docs] def getFieldNames(self):
"""
Gets the field names for the data class
:rtype: list
:return: list of attribute names for the data class
"""
return [x.name for x in fields(self)]
[docs] def getValues(self):
"""
Gets the values for attributes that are not none
:rtype: dict
:return: dictionary where key is the attribute name value is the value
of the attribute
"""
values_dict = {
x: getattr(self, x)
for x in self.getFieldNames()
if getattr(self, x) is not None
}
return values_dict
[docs]class HelixArcAtomGroup(MSAtomGroup):
"""
Class containing group of HelixArc and combine then into helix groups
"""
RADIUS_TOLERANCE = 20
BEND_TOLERANCE = 40
BASE_ATOM_TYPE = HelixArc
[docs] def __init__(self, struct, group_atom_ids):
"""
Initialize HelixArcAtomGroup
:type struct: `schrodinger.structure.Structure` or `cms.Cms`
:param struct: Structure to which the atom group belongs. It will be
cms model if trajectory is used, and structure if trajectory is not
used.
:type group_atom_ids: list
:param group_atom_ids: list of atom ids to add to the group
"""
super().__init__(struct, group_atom_ids, GROUP_CENTROID)
[docs] def estimateRadius(self, arcs, all_pos):
"""
Estimate helix radius from arc points by finding their distance from
the helix axis.
:param arcs: The arcs making the helix
:type arcs: list of HelixArc
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
:returns: mean radius for the helix
:rtype: float
"""
# Get atom positions in the helix, avoiding repeating positions to
# remove bias
helix_pos = []
for arc in arcs:
arc_pos = [tuple(x) for x in arc._getAllXYZ(all_pos).tolist()]
helix_pos.extend(arc_pos)
helix_pos = np.array(msutils.get_unique_ordered_list(helix_pos))
# Center at origin
helix_center = helix_pos.mean(axis=0)
helix_pos -= helix_center
# Fit 3D line to helix points, the first vector of vh is the direction
# of the helix
_, _, vh_val = np.linalg.svd(helix_pos)
helix_axis = vh_val[0]
# Calculate perpendicular distance of each point from the helix axis
radii = []
for pos in helix_pos:
v_perpendicular = pos - np.dot(pos, helix_axis) * helix_axis
radii.append(transform.get_vector_magnitude(v_perpendicular))
return np.mean(radii)
[docs] def fixDiameter(self, helix_grp, all_pos):
"""
Estimate the diameter of helix group and then set the helix group
diameter as the diameter of each arc. Since the local diameter of arc
may not truly represent the helix diameter when the helix is not
perfectly straight or points are far apart.
:param helix_grp: The list of arcs making the helix group
:type helix_grp: list
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
"""
narcs = len(helix_grp)
base_per_turn = int(np.mean([arc.base_per_turn for arc in helix_grp]))
# Split the helix groups in turns to find radius of each turn
# This allows for bending of helix
if base_per_turn < narcs:
split_arcs = [
helix_grp[x:x + base_per_turn]
for x in range(narcs - base_per_turn + 1)
]
else:
split_arcs = [helix_grp]
# Hold arc diameter from each split
for arc in helix_grp:
arc.diameter = []
# Estimate the radius of each subgroup and set that as the local arc
# diameter. We do not record radius of a helix since diameter is the
# convention
for arcs in split_arcs:
radius = self.estimateRadius(arcs, all_pos)
for arc in arcs:
# Add diameter from each split and average later
arc.diameter.append(2 * radius)
# Average arc diameter from each split
for arc in helix_grp:
arc.diameter = np.mean(arc.diameter)
[docs] def groupAtoms(self):
"""
Overwrite parent group to create group where each group element is
a HelixArc comprising for four atoms
"""
self.group_atoms = []
num_atoms = len(self.group_atom_ids)
if num_atoms < 4:
raise AttributeError('ArcGroup can not be created for less '
'than 4 atoms.')
for index in range(num_atoms - 3):
atom_indices = self.group_atom_ids[index:(index + 4)]
self.group_atoms.append(
self.BASE_ATOM_TYPE(atom_indices, self.struct, GROUP_CENTROID))
[docs] def appendToSplitGroup(self, split_grp, current_grp):
"""
Append the current group to split group if current group contains at
least three atom group(s).
:type split_grp: list
:param split_grp: list of list of atom group(s)
:type current_grp: list
:param current_grp: list of atom group(s) to append
"""
if len(current_grp) > 3:
split_grp.append(current_grp)
[docs] def cleanBoolCheck(self, bool_check):
"""
Clean the check to remove a single noise False point
:type bool_check: list
:param bool_check: list of bool values
"""
for index in range(1, len(bool_check) - 1):
if bool_check[index]:
continue
if bool_check[index - 1] and bool_check[index + 1]:
bool_check[index] = True
return bool_check
[docs] def splitGroupByChecks(self, groups, grp_check):
"""
Split an ordered list of elements into list of list containing the
(not all) elements of arc_group. Splitting occurs when False values are
encountered in grp_check, and the corresponding groups elements are
excluded.
:type groups: list
:param groups: list of elements to split
:type grp_check: list
:param grp_check: list of bool values
:rtype: list(list)
:return: list containing list of elements whose corresponding index in
grp_check is True
:raises AttributeError: If the length of groups and grp_check is not
the same
"""
if len(groups) != len(grp_check):
raise AttributeError(
'Length of group list is not same as as the checking list.')
# Split the groups into split groups
split_groups = []
current_group = []
for grp_element, allowed in zip(groups, grp_check):
if allowed:
current_group.append(grp_element)
else:
self.appendToSplitGroup(split_groups, current_group)
current_group = []
self.appendToSplitGroup(split_groups, current_group)
return split_groups
[docs] def splitGroupByArc(self, groups):
"""
Split the groups into multiple list of groups that contain atom group
that are arcs
:type groups: list
:param groups: list of atom group
:rtype: list(list)
:return: list of list of atom groups that are arc
"""
is_arc = [x.is_arc for x in groups]
return self.splitGroupByChecks(groups, is_arc)
[docs] def splitGroupByLine(self, arc_grps):
"""
Split arc groups into helix groups where arc groups should lie in a
straight line. The amount of straight-ness is controlled by bend
tolerance.
:type arc_grps: list
:param arc_grps: list of atom groups that are arc
:rtype: list
:return: list of atom groups that are helical arcs
"""
split_groups = []
for arc_grp in arc_grps:
bend_check = []
for index in range(1, len(arc_grp)):
# Get vector with last arc center
vector = arc_grp[index].center - arc_grp[index - 1].center
arc_grp[index].vector = transform.get_normalized_vector(vector)
if index == 1:
# Assuming first two arcs to be in line
bend_check = [True] * 2
else:
# Check vector straight-ness
dot_val = np.dot(arc_grp[index].vector,
arc_grp[index - 1].vector)
bend_angle = np.rad2deg(np.arccos(abs(dot_val)))
bend_check.append(bend_angle < self.BEND_TOLERANCE)
bend_check = self.cleanBoolCheck(bend_check)
split_groups.extend(self.splitGroupByChecks(arc_grp, bend_check))
return split_groups
[docs] def splitGroupByRadius(self, helix_grps, all_pos):
"""
Split helix arc groups into helix arc groups where local arc radius are
similar. The amount of similarity is controlled by radius tolerance.
:type helix_grps: list
:param helix_grps: list of atom groups that are arc
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
:rtype: list
:return: list of atom groups that are helical arcs
"""
# Find radius of helix from current groups
for helix_grp in helix_grps:
self.fixDiameter(helix_grp, all_pos)
split_groups = []
for helix_grp in helix_grps:
grp_radius = [x.diameter / 2 for x in helix_grp]
avg_radius = trim_mean(grp_radius, self.RADIUS_TOLERANCE / 100)
radius_check = [
abs(x - avg_radius) * 100 / avg_radius < self.RADIUS_TOLERANCE
for x in grp_radius
]
radius_check = self.cleanBoolCheck(radius_check)
split_groups.extend(self.splitGroupByChecks(helix_grp,
radius_check))
return split_groups
[docs] def helixGroupToData(self, helix_grp):
"""
Get the average helix property from list of local helix arcs
:type helix_grp: list
:param helix_grp: list of local helix arc
:rtype: `HelixData`
:return: Helix data containing helix properties
"""
# Set local properties from individual arcs
values = {}
for prop in HelixData().getFieldNames():
prop_vals = [getattr(x, prop, None) for x in helix_grp]
if any(prop_vals):
values[prop] = np.mean(prop_vals)
# Create string of aids
aids = set()
for arc in helix_grp:
aids.update(arc.indices)
aids_str = HELIX_AIDS_DELIMITER.join(map(str, aids))
return HelixData(**values, narcs=len(helix_grp), aids=aids_str)
[docs] def calculateHelices(self, all_pos=None):
for arc in self.group_atoms:
arc.calculateHelixProp(all_pos)
# Group the arcs
helix_grps = self.splitGroupByArc(self.group_atoms)
# Group the arcs by in line
helix_grps = self.splitGroupByLine(helix_grps)
# Group the arcs by radius
helix_grps = self.splitGroupByRadius(helix_grps, all_pos)
# Convert to dataclass
helix_data = [self.helixGroupToData(x) for x in helix_grps]
return helix_data
###############################################################################
# Atom Groups - Constructed using multiple MSAtomGroup
###############################################################################
[docs]class MSAtomGroupsBase:
"""
Class containing atom groups and calculate properties for each of them
"""
BASE_ATOM_GROUP = MSAtomGroup
[docs] def __init__(self,
struct,
atom_ids=None,
msys_model=None,
frames=None,
wrap_traj=False,
logger=None):
"""
Initialize MSAtomGroups
:type struct: `schrodinger.structure.Structure` or `cms.Cms`
:param struct: Structure to which the atom group belongs. It will be
cms model if trajectory is used, and structure if trajectory is not
used.
:type atom_ids: list or None
:param atom_ids: atom indices to create atom groups from
:type msys_model: `msys.System` or None
:param msys_model: Desmond model system
:type frames: list or None
:param frames: list of trajectory frames if analyzing trajectory. None
if analyzing static structure.
:type wrap_traj: bool
:param wrap_traj: If true, it will wrap trajectory
:type logger: function or None
:param logger: logger function to log updates to.
"""
self.struct = struct
self.msys_model = msys_model
self.frames = frames
self.wrap_traj = wrap_traj
self.logger = logger
if atom_ids:
self.atom_ids = atom_ids
else:
self.atom_ids = list(range(1, struct.atom_total + 1))
self.groups = []
self.setAtomGroups()
[docs] def updateFrameLog(self, frame_id, prop_name):
"""
Add to log with updates to calculations
:param int frame_id: Index of the current frame id
:param str prop_name: Name of the property being calculated
"""
if not self.logger or not prop_name:
return
msg = (f'Calculating {prop_name} for frame {frame_id} of '
f'{len(self.frames)}...')
self.logger(msg)
[docs] def addIndicesToGroup(self, atom_indices):
"""
Create a group from passed atom indices and add it to groups
:type atom_indices: list
:param atom_indices: list of atom indices
"""
self.groups.append(self.BASE_ATOM_GROUP(self.struct, atom_indices))
[docs] def setAtomGroups(self):
"""
Function to create atom groups. Runs when the object is initiated.
:raises NotImplementedError: Raises error if not implemented for the
class
"""
raise NotImplementedError(
'Atom group creator has not been implemented.')
def _getCords(self, pos=None):
"""
Get xyz values for the all the atom groups
:type pos: `numpy.ndarray` or None
:param pos: all coordinates in the trajectory for the current/analyzed
frame. None if it is a static structure
:rtype: list
:return: xyz values for all the atom groups
"""
return [x.getGroupXYZ(pos) for x in self.groups]
[docs] def getCordsIterator(self, prop_name=None):
"""
Iterator for xyz for all the groups. This will iterate over all
trajectory frames or once for a static structure
:type prop_name: str or None
:param prop_name: Name of the property which is being calculated. If
logger is available and property is not None, it will log the
progress.
:rtype: list
:return: list of xyz values for all groups
"""
if self.frames:
for frame_id, frame in enumerate(self.frames, 1):
self.updateFrameLog(frame_id, prop_name)
if self.wrap_traj:
pbc = analysis.Pbc(frame.box)
coordinates = frame.pos()
# Wrap the trajectory coordinates to first box.
coordinates[:] = pbc.calcMinimumImage(
np.zeros([1, 3]), coordinates)
else:
topo.make_whole(self.msys_model, [frame])
coordinates = frame.pos()
yield self._getCords(coordinates)
else:
yield self._getCords()
[docs]class MSMoleculeAtomGroups(MSAtomGroupsBase):
"""
Class containing atom groups created for each molecule in the structure
"""
[docs] def setAtomGroups(self):
"""
Set atom groups for each segment in the structure
"""
for mol in self.struct.molecule:
atom_indices = [
x for x in mol.getAtomIndices() if x in self.atom_ids
]
if atom_indices:
self.addIndicesToGroup(atom_indices)
[docs]class HelixGroups(MSAtomGroupsBase):
"""
A class containing group of HelixArcAtomGroup and calculate global helix
properties from it
"""
BASE_ATOM_GROUP = HelixArcAtomGroup
[docs] def setAtomGroups(self):
"""
Set atom groups for each segment in the structure
"""
st_graph = analyze.create_nx_graph(self.struct)
for segment in graph.find_segments(st_graph):
atom_indices = [x for x in segment if x in self.atom_ids]
# Need at least 2 full dihedral in the assumption
if len(atom_indices) < 8:
continue
self.addIndicesToGroup(atom_indices)
[docs] def getHelixData(self, all_pos=None):
"""
Calculate all helices for current frame / structure
:type all_pos: `numpy.ndarray` or None
:param all_pos: all coordinates in the trajectory for the
current/analyzed frame. None if it is a static structure
:rtype: list
:return: list of HelixData that contain properties of each helix
"""
all_helix_data = []
for group in self.groups:
all_helix_data.extend(group.calculateHelices(all_pos))
return all_helix_data
[docs] def getHelixDataIterator(self):
"""
Iterator to get helices in each frame. Will yield a single list of
helices for a static structure.
:rtype: list
:return: list of all the helix data containing helix properties
"""
if self.frames:
for frame_id, frame in enumerate(self.frames, 1):
topo.make_whole(self.msys_model, [frame])
coordinates = frame.pos()
all_helix_data = self.getHelixData(coordinates)
# Update data with current time
for helix_data in all_helix_data:
helix_data.time = frame.time / 1000
yield all_helix_data
else:
yield self.getHelixData()