symmetry module

Expansion of molecular & crystal symmetry.

Direct stand-alone usage has been deprecated; see symexp.py.

Module author: Michael S. Chapman <chapmanms@missouri.edu>

Authors:

Michael S. Chapman <chapmanms@missouri.edu>

Oregon Health & Science University & University of Missouri

Changed in version Start: 09/09/11

Changed in version 0.3.7: 11/04/11

Changed in version 0.5.0: 05/04/15 ReStructuredText docs

Changed in version 1.0.0: 11/06/20 Python 2.7 –> 3.6

Changed in version 1.0.8: 01/20/25 catch roundoff error for acos(), line 206

Changed in version 1.1.1: 10/17/25 complete supports addition of residues to neighbor oligopeptides.

class symmetry.Arguments(imports=[], main=None, *args, **kwargs)

Bases: Arguments

Manager for command-line arguments from main and imported modules.

This is a template to be copied into modules, for the handling of command- line options.

methods export() and domestic() should be overridden by the module subclass with declarations of the command-line arguments needed by the module.

#QQQ Should arguments be forced to those compatible with energy?

Parameters:

  • imports ((list of) ArgumentParser(s)) – (list of) module(s) from which to obtain “parent” ArgumentParser objects for inclusion.

  • main (bool|NoneType) – Add information and arguments appropriate for a main program (or not); if None, will determine by whether this class is defined in __main__.

domestic()

Defines options used only when main program and not when imported.

export()

Defines options used in both stand-alone and imported modes.

class symmetry.Cartesian_operator(name='', rotation=[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], translation=[0.0, 0.0, 0.0], verbose=True)

Bases: object

Rotation + translation in Cartesian space.

x’ = [R]x + t.

Parameters:

  • rotation (numpy.matrix-like.shape(3,3)) – R - row-major 3x3 rotation matrix.

  • translation (numpy.matrix-like.shape(3,1)) – t - translation vector applied after rotation.

  • name (str) – optional name, 1st character used chain ID, 1st 4 characters used for PDB segid.

  • verbose (bool) – print information about operator.

euler(verbose=False)

Euler rotation and axis from rotation matrix.

Parameters:

verbose (bool) – print rotation & axis

Returns:

(rotation (degrees), unit-vector-axis)

Return type:

float, numpy.matrix.shape(3,1)

isunit()
Returns:

True if unit operator.

Return type:

bool

class symmetry.Documentation

Bases: object

User-level documentation (imports).

class Input

Bases: object

property local

Generic documentation for import into other modules.

class symmetry.Symmetry(local=None, unit=True, space_group=None, cell=None, lattice_translations=1, distance=3.5, tolerance=0.5, completion=None, capping=False, verbose=False)

Bases: object

Symmetry expansion of coordinates.

Changed in version 0.5.5: 10/20/15 Riding anisotropic Us supported. Tested by comparing traces of ANISOU PDB records for symmetry equivalents, and visual inspection of R3 expanded insulin.

For local symmetry, roughly O(NlogN) to find neighbors:

  • 25 ms for 50 atom target

  • 15 s for 4200 atom target

Re-expansion O(1), ~1 ms.

Variables:

  • operators (tuple(int,int,int,int,int)) – indices of the operators used to generate each neighbor from unique (local, crystallographic, translation_a, b, c), noting that only local currently implemented.

  • ~.neighbors (list(ndarray(type=bool))) – selection of atoms from unique that generate neighbors with the corresponding transformation from operators.

Todo

Merging equivalents back into protomer. Symmetry has been designed on the premise that there is a complete protomer from which equivalents are generated by constrained operators. A desirable extension would be to merge changes in equivalents back into the protomer, supporting refinement of volumes that might contain symmetry-equivalent regions of interest, but not fully enclosing one whole protomer. There will need to be user-driven algorithms for precedence or weighting of equivalents when more than one is changing. These should probably include a user-designated favorite protomer (could be multi-chain); weighting by distance from center of mass; or by bond count to fixed neighbors.

Initialize operators for expansion / set default neighbor criteria.

Parameters:

  • local (list of N (Carterian_operator or tuple for instantiation of Cartesian_operator).) – local (molecular) operators or None.

  • unit (bool) – insert (if absent) and use the unit operator.

  • space_group (str or NoneType) – (optionally shortened) Hermann-Mauguin symbol with ‘:setting’ appended (if required; see space_group, or None.

  • cell (tuple(6) or numpy.matrix-like (size = 3x3)) – unit cell parameters, or fractional to orthogonal (Å) matrix, required if space_group.

  • lattice_translations (int) – Number to search in each direction for neighbors. If unique atoms are within -0.5 to +1.5 (fractional units, most cases) then 1 should suffice. In some cases a larger search generates more neighbors - wise to check 2. Zero disables.

  • distance (float) – within which neighboring atoms will be found, see also Symmetry.neighbors.

  • tolerance (float) – neighbors to be re-identified when an atom has moved more than this amount (Å), see also Symmetry.neighbors.

  • completion (string or NoneType) – if ‘residue’ or ‘chain’, expand the selection of neighboring atoms to the whole group.

  • capping (bool) – if completion = ‘residue’, for gaps in equivalents substitute NME and ACE caps for the preceding and following residue so that artificial charges are neutralized (for force field restraints).

  • verbose (bool) – print diagnostics.

Warning

This instantiation assumes that the non-crystallographic symmetry, unit cell and magnification will remain constant. Refinement will require re-instantiation. RSRef currently assumes that magnification will only be refined for EM which where space_group=None, so that it is safe to instantiate only at start up.

Note

local could be the result of symmetry.safe_eval of an input tuple of which the following is a short example:

("p",        ((0.500000,       0.809000,      -0.309000),
             (-0.809000,       0.309000,      -0.500000),
             (-0.309000,       0.500000,       0.809000)),
              (0.000000,       0.000000,       0.000000)),
("e",        ((0.309000,      -0.500000,       0.809000),
              (0.500000,       0.809000,       0.309000),
             (-0.809000,       0.309000,       0.500000)),
              (0.000000,       0.000000,       0.000000)))
class Cap(atoms, neighbors, verbose=True)

Bases: object

Conversion of 1st & last residues of neighbor fragments to NME, ACE caps

When using force field (stereochemical) restraints, N-terminal and C-terminal ends are charged (NH3+ and COO- respectively). If these termini are at the ends of neighbor fragments generated from symmetry, then electrostatic artifacts may be introduced. These can be mititaged by neutralizing NH3+ with an acetyle group (ACE) and COO- with an N-methylacetamide group (CH3-CO-NH-CH3), better mimicking the environment within the integral protein.

Be careful, because some Google searches incorrectly assert that NME caps the N-terminus and ACE the C-terminus.

In our implementation, ACE and NME are substituted for the 1st and last residues of a neighbor fragment on the assumption that fragments have already been extended by >= 1 residue beyond those of user-interest. This then allows cap atom positions to be updated according to changes in the positions of most corresponding protomer atoms. This algorithm is applicable only to symmetry equivalents and not to the original protomer, but this is moot, because the program does not know, for the protomer, which terminii are natural and which represent the truncated ends of what is experimentally resolvable. The philosophy here is that the user should cap unnatural truncations in the input PDB file and that support is only needed for the capping of symmetry-generated fragments.

Atom name substitutions can, at best, approximate the appropriate stereochemistry. For example, there are C–N bonds being changed to C–H without change to the bond length. It is not clear yet whether this approximation will be adequate or whether regularization will be needed.

It is also not yet known whether PDBfixer will require either all or no hydrogens or whatever is available. Also not known whether a TER record will be needed before ACE.

ACE & NME residue names might be specific for Amber force fields.

Variables:

  • aa2ace – Atom name translations to convert amino acid atom names to those of corresponding cap atoms. Attention has not (yet) been paid to implied chirality, if this is important - it might require swapping of some of the equivalences.

  • aa2nme – Atom name translations to convert amino acid atom names to those of corresponding cap atoms. Attention has not (yet) been paid to implied chirality, if this is important - it might require swapping of some of the equivalences.

For a fragmented symmetry equivalent, substitute caps for terminal residues

Actually, figure out how on calling method cap.

Parameters:

  • atoms (Atoms) – Protomer coordinate set with chain, (segid) and residue names/numbers that correspond to neighbors. This will usually be the untransformed protomer.

  • neighbors (ndarray-like, dtype=bool, size=# atoms.) – True for atoms that are neighbors in a symmetry equivalent.

Variables:

  • ncap – list of atoms indices for the CA of each residue to be N-terminally capped where neighbors indicates fragmented equivalents.

  • ccap – list of atoms indices for the CA of each residue to be C-terminally capped where neighbors indicates fragmented equivalents.

aa2ace = {'C': 'C', 'CA': 'CH3', 'CB': 'H2', 'HA': 'H1', 'HA2': 'H1', 'HA3': 'H2', 'N': 'H3', 'O': 'O'}
aa2nme = {'C': 'H3', 'CA': 'C', 'CB': 'H2', 'CD': 'H', 'H': 'H', 'HA': 'H1', 'HA2': 'H1', 'HA3': 'H2', 'N': 'N'}
mutate(equivalent, neighbors, verbose=True)

Change the residue name & atom attributes to convert aa to cap.

Parameters:

  • equivalent (Atoms) – A copy of the unique protomer, before or after it has been operator-transformed.

  • neighbors (ndarray-like, dtype=bool, size=# atoms.) – True for atoms that are neighbors in equivalent after operator-transformation. Atoms unused after conversion to the cap should have already been set to False using self.prune().

Returns:

a mutated deep copy of equivalent in which the selected atoms of relevant amino acids have residue & atom names changed to the cap forms. Amino acid atoms that will not be used in the cap remain with unchanged names, but should be marked for non- use by self.prune().

Note that a Cap object is specific for a symmetry equivalent

Warning

Supports only Amber-like atom naming convention.

prune(neighbors, verbose=True)

Remove amino acid atoms that do not correspond to cap atoms.

Parameters:

neighbors (ndarray-like, dtype=bool, size=# atoms.) – True for atoms that are neighbors in a symmetry equivalent.

Warning

Supports only Amber-like atom naming convention.

complete(atoms, neighbors, completion=None, extend=1)

Expand list of neighbors to include all atoms w/in residue or chain.

Requires atoms in the same residue to share the same chain and segids (even if null) and for at least one of these to be unique for the residue number.

Parameters:

  • atoms (Atoms) – Protomer coordinate set with chain, (segid) and residue names/numbers that correspond to neighbors. This will usually be the untransformed protomer.

  • neighbors (ndarray-like, dtype=bool, size=# atoms.) – True for atoms that are neighbors in a symmetry equivalent.

  • completion (str or NoneType) – ‘chains’, ‘residues’, ‘atoms’ aka None.

  • extend (int) – number of residues to add, each end, to residue-completed neighbors. This should be 1 (or more) for OpenMM use, to avoid exception for single amino acid neighbor fragments.

Returns:

completed neighbors.

Return type:

ndarray-like, dtype=bool, size=# atoms.

Warning

Compatibility of completion options: completion=None will be incompatible with energy restraints; completion=’residue’ will leave symmetry-related oligopeptides with charged N- and C- terminal ends that can distort energy restraints. For example, a starting structure for the ~500 amino acid AAV capsid protein, has an energy of 36,000 kJ/mol isolated, an approximately proportional 400,000 kJ/mol as 10 complete chains, but 1,450,000 kJ/mol as a single subunit surrounded by fragments equivalent, in total size, to about 1.5 additional subunits. The additional 1M kJ/mol coming from the artifactual charges are likely to spoil refinement, even if these atoms are fixed. Capping will be needed to avoid this.

Todo

Capping: Need to add the capability for ACE capping of “unnatural” amino termini and NME capping of carboxylate fragment termini. Anticipating that the best approach will be to use complete() to extend one or more residues, then replace the new ending residues with ACE or NME. Note that the capping atoms will not have direct correspondence with those of the unique protomer.

debug = False
equivalent(xyz, ilocal=0, ilattice=0, a_trans=0, b_trans=0, c_trans=0, anisou=None)

Transform coordinates by a single set of local and lattice operators.

Parameters:

  • xyz (ndarray.shape(3,N)) – coordinates (Angstrom).

  • ilocal (int) – index of Symmetry.local operator to be applied.

  • ilattice (int) – index of Symmetry.lattice operator to be applied (not implemented).

  • a_trans (int) – unit cell translations along a axis.

  • b_trans (int) – unit cell translations along a axis.

  • c_trans (int) – unit cell translations along a axis.

  • anisou (Anisotropics/NoneType) – U tensors if available, else None.

Returns:

transformed xyz, anisou

Return type:

ndarray, Anisotropics

inverse(i)

Identify the inverse of a local operator, if present.

Parameters:

i (int) – index of original local symmetry operator

Returns:

index of inverted operator

Return type:

int|NoneType

neighbors(unique, target_atoms=None, selection=None, distance=None, lazy=True, tolerance=None, completion=None)

Symmetry-expand unique for all neighbors of selected atoms in target.

Parameters:

  • unique (Atoms) – atomic coordinates to be symmetry expanded.

  • target_atoms (Atoms or NoneType) – atoms whose neighbors from unique (and unique’s symmetry equivalents) will be found. If None, defaults to unique.

  • selection (ndarray, dtype=bool) – subset of target atoms to use. True to select all. Defaults to target_atoms.evaluate if target_atoms defined, or True.

  • distance (float) – limit in Å; if None, use instantiated value.

  • lazy (bool) – recalculate which atoms are neighbors only if the selection unique has changed by more than tolerance

  • toleranceprune (float) – see lazy; if None, use instantiated value.

  • completion (str or NoneType) – if ‘residue’ or ‘chain’, expand the selection of neighboring atoms to the whole group.

Returns:

(input atoms w/ unique.close reset, symmetry-related neighbors) unique.close will be True for atoms w/in distance of selection and not unique.evaluate. unique retains all atoms and is suitable for PDB output. Arrays of the symmetry equivalents are compressed, containing only neighboring atoms, and not necessarily whole residues or molecules suitable for output etc.

Return type:

(Atoms, Atoms)

old_with_neighbors(unique)

Iterate over selected atoms and (symmetry-related) neighbors.

Equivalent to Symmetry.neigbors (with default arguments), then iterating over (Atoms.evaluate | Atoms.close) for asymmetric unit and equivalent atoms.

Parameters:

unique (Atoms) – atomic coordinates, (NCS) asymmetric unit.

Returns:

atom type, position, occupancy, b-factor, anisotropic-Us, evaluate, close, atom index (from which equivalent generated), operators=(local, lattice, trans_a, trans_b, trans_c).

Return type:

string, NumPy.matrix(dtype=float, shape=(3,1)), float, float, NumPy.matrix(dtype=float, shape=(3,3))/NoneType, bool, bool, int, tuple(5 int)

Changed in version 1.1.0: (10/6/25) yields the atoms atom-index instead of the atom number with which it was associated in the input PDB.

rechain(expanded, fragment=False)

Rename chains by symmetry operator / fragment

Parameters:

  • expanded (Atoms)

  • fragments (bool) – New chain ID where there is a break in polymer residue number (needed for OpenMM systems).

Returns:

expanded with updated chain IDs

Return type:

Atoms

The kernel is Atoms.ChainName which has not been retested thoroughly, now adding symmetry operator to chain for the attributes that should trip generation of a new chain ID.

rechain gives unique chain identifiers to each (peptide) fragment of a symmetry equivalent, thereby ensuring that openmm forcefields do not connect them with a covalent bond. This should not be necessary if the ends of fragments are capped or have terminus-specific residue names (NALA, CALA, NARG, CARG…). rechain() should be mostly benign, but should not be used unless necessary: (1) Chain name substitution is confusing to users, particularly when also combined with openmm renumbering of residues from 1 with every chain. (2) Possibility that the full range of single character IDs will be exceeded. This method will use lower case, numeric and punctuation, but there is no guarantee that the extended set will be supported by external modules.

repr_op(op_indices)

Self-documenting representation of the operator set.

Parameters:

op_indices (tuple(int,)*5) – indices of the local & lattice operators and unit translations along each axis.

verbose = 1
with_neighbors(unique)

Iterate over selected atoms and (symmetry-related) neighbors.

Equivalent to Symmetry.neigbors (with default arguments), then iterating over (Atoms.evaluate | Atoms.close) for asymmetric unit and equivalent atoms.

Parameters:

unique (Atoms) – atomic coordinates, (NCS) asymmetric unit.

Returns:

atom type, position, occupancy, b-factor, anisotropic-Us, evaluate, close, atnam, altloc, resnam, chain, segid, resnum, insert, element, charge, polymer operators=(local, lattice, trans_a, trans_b, trans_c).

Return type:

string, NumPy.matrix(dtype=float, shape=(3,1)), float, float, NumPy.matrix(dtype=float, shape=(3,3))/NoneType, bool, bool, string, string, string, string, string, int, string, string, string, tuple(5 int)

Changed in version 1.1.0: (10/6/25) yields the atoms atom-index instead of the atom number with which it was associated in the input PDB.

Changed in version 1.1.1: (10/28/25) Replace output of atom-index with explicit atnam, altloc, resnam, chain, segid, resnum, insert, element & charge. (It has become too difficult to maintain (more efficient?) referal to the unique protomer for these attributes with introduction of methods (capping) that insert and delete atom records.)

class symmetry.Target(xyz, use=None, distance=3.5)

Bases: object

Atomic coordinates for which neighbors are sought.

Sphere and box that enclose target atoms and anything w/in distance.

Parameters:

  • xyz (Atoms) – atomic coordinates whose neighbors are to be found.

  • use (ndarray, dtype=bool) – subset of atoms to use, defaults to all.

  • distance (float) – limit in Angstrom for neighbors to be considered close.

debug = False
isclose(xyz, distance_increase=0.0)

Is point w/in distance of target sphere?

Parameters:

  • xyz (numpy.matrix(shape(3,1), dtype=float)) – position, Cartesian Angstrom coordinates.

  • distance_increase (float) – to be added to self.distance.

Returns:

close

Return type:

bool

isneighbor(xyz)

Which of N points in xyz is w/in self.distance of any self (target) atom?

Note that this does not care about atomic properties at xyz, including whether an atom has been selected for evaluation. However, the target coordinates (for which neighbors are sought) will have previously been subject to selection criteria.

Parameters:

xyz (numpy.matrix(shape(3,N), dtype=float)) – position, Cartesian Angstrom coordinates.

Returns:

close

Return type:

ndarray(shape=N, dtype=bool)

symmetry.safe_eval(filename)