restraints module

Ancilliary restraints (B-factor).

Module author: Michael S. Chapman <chapmami@ohsu.edu>

Module author: Leo Selker

Authors:

Michael S. Chapman <chapmami@ohsu.edu>,

Oregon Health & Science University

Version:

0.5, March 23, 2016

Changed in version 0.1.1: 08/29/11 Start

Changed in version 0.5.0: (5/2/15) ReStructured Text documentation

Changed in version 0.5.6: (10/21/15) Adding vdW restraint from Leo

See also

Limit() in torsion.py.

class restraints.Arguments(imports=[], main=False, *args, **kwargs)

Bases: argparser.ArgumentParser

Command-line arguments for optimization restraints

Parameters:

• imports (list) – method(s) adding arguments to ArgumentParser instance. (Alternative to ArgumentParser parents w/ better maintained groups.)
• main (bool) – Parser for main program (i.e. not listed within an ArgumentParser parents argument). Adds default program information.

static export()

class restraints.Brestraint(atoms, weight, group_exclusions=None)

Bases: object

Restraint on variation of adjacent atomic displacement parameters.

Parameters:

• atoms (Atoms) – atomic parameters (coordinates etc).
• weight (float) – empirical (user-adjusted) weight for the restraint.
• group_exclusions (class Group, bool or NoneType) – intra-selection restraints are excluded within this Group of (rigid) Selections.

df(b)

Atomic partial derivatives for B-factor variation restraint.

Parameters:b (ndarray(type = float)) – current atomic displacement parameters (must be aligned w/ self._bonded) Returns:weight * ∂Σ_j [B_i - B_j]²/∂B_j, for all j neighbors of i atoms Return type:ndarray(type = float)

f(b)

Objective function for B-factor variation restraint.

Parameters:b (ndarray(type = float)) – current atomic displacement parameters (b must be aligned with self._bonded) Returns:Σ_{i, j} [B_i - B_j]² * weight, for all j neighbors of i atoms Return type:float

rms(b)

RMS variation in B-factor for bonded atoms.

Parameters:b (ndarray(type = float)) – current atomic displacement parameters (b must be aligned with self._bonded) Returns:(Σ_{i, j} [B_i - B_j]²)^0.5, for all j neighbors of i atoms Return type:float

class restraints.VDWRestraint(atoms, weight=1, parameterization=None, symmetry=None, max_shift=1.0, radiiPath=None, exclude_ss=True, same_conformer=True, conformer_coupling=0.0)

Bases: object

Restraint on overlapping van der Waal radii of atoms

Excludes all pairs of SG atoms on basis that might be disulfide.

Changed in version Reimplemented: to reduce incapacitating memory usage. Previously stored several (natom, natom) float matrices, and now stores only boolean matrices designating which interactions to recalculate, now every time. Nearly 100x faster.

Warning

Does not support “open” symmetry, only closed symmetry in which repeated application of the same operator returns you to the starting point. Algorithms use efficiencies of closed symmetry in searching for potential contacts. Closed symmetry includes all crystal symmetry and most molecular symmetry, but not all non-crystallographic symmetry.

Warning

applies the full restraint to atoms with partial occupancies, ignoring the possibility that in multi-conformer structures, conformers might be selected in mutually compatible (non-overlapping) ways.

Todo

account for open symmetry

Todo

account for states that are only sometimes occupied

Parameters:

• atoms (Atoms) – coordinates (used to extract atom types etc., but positions will be updated for each restraint calculation).
• weight (float) – multiplier for the vdW residual restraint.
• parameterization (:py:class:’AtomicParameterization’ or sub-class) – definition of what parameters are refining and how they depend on atomic parameters. This is used here to limit calculation of vdW overlap & partials to those that will affect moving atoms.
• symmetry (Symmetry/NoneType) – defines molecular & lattice symmetry, if any.
• radiiPath (str) – file name for a python dictionary to extend or over-ride default vdW radii.
• max_shift (float) – maximum expected change in atomic coordinates, used in pre-calculating possible contacts.
• exclude_ss (bool) – exclude contacts between pairs of SG atoms on the assumption that they are CYS residues and, if close, they are likely disulfide-bridged.
• same_conformer (bool) – overlap is only considered between atoms from the same conformer (A, B, ...) or if one of the pair is single- conformer (no letter).
• conformer_coupling (float) – NOT IMPLEMENTED YET... if 1.0, minimal overlap will be calculated, assuming that alternate conformers (where present) completely resolve any overlap. This could be a reasonable approximation if “common atom” models and pseudo-symmetry in rotamers were considered (neither currently supported). If conformer_coupling=0.0, conformer selections are assumed to be uncoupled (random) and therefore the restraint is multiplied by the product of the two atoms’ occupancies. Both are approximations, coupled likely to underestimate vdW overlap and uncoupled likely to overestimate. Intermediate values 0.0 < conformer_coupling < 1.0 will lead to a weighted sum of the two approaches. Note that inappropriate (un)coupling could lead to over- or under- estimation of the conformational diversity if their occupancies are being refined.

Variables:

distance – cut-off used to determine from the initial structure which atom pairs might get close enough for van der Waals repulsion. This distance is also used as the cut-off for the generation of symmetry equivalents, unless the equivalents have already been generated with a larger cut-off.

..warning:: If the degree of coupling is not appropriately approximated,

then refinement of conformer occupancies can be biased by over- or under-estimating the degree of real van der Waals overlap. Should consider whether really need to use the restraint when refining conformer occupancies.

..warning:: Currently considering refinement of the occupancies of

rotamers generated on the fly, but this class will only consider positions explicitly defined in the atomic coordinates.

Todo

Exclude contacts from different PDB-defined conformers.

Todo

Exclude contacts w/in same residue, assuming that they will be handled by stereochemical constraints / restraints.

Todo

Handling of hydrogens. The default radii assume that hydrogens are explicitly defined, but most (low resolution) structures do not include them.

class Status(parent, atoms)

Bases: object

Status of vdW contacts.

Parameters:

• parent (VDWRestraint) – self of the owning class
• atoms (Atoms) – current coordinates

VDWRestraint.connectivityExclusion = 3

Atoms connected by this or fewer bonds will be excluded from the restraint.

E.g. if atoms a, b, c, d are bonded: a–b–c–d, then all pair except a ... d will be excluded.

In other’s refinements, torsion angles would be restrained explicitly, not through vdW (connectivityExclusion = 3), and even without torsional restraints a value of 3 might be needed to avoid locking dihedrals to exactly staggered conformations. With a value of 2, expect dihedrals to be regularized.

VDWRestraint.distance = 5

Neighbors within at least distance will used in vdW calculations.

The actual distance might be larger, if symmetry-equivalent neighbors have previously been calculated with a larger distance.

VDWRestraint.f(atoms)

Objective function for VDW restraint

Parameters:atoms (Atoms) – current (trial) coordinates Returns:V = w Σ:sub:’i|isin|N’ Σ:sub:’j|isin|N’ {(|**x**:sub:’i’ - |x:sub:’j’| - (r:sub:’i’ + r:sub:’j’)}:sup:’p’, where V is the van der Waals penalty summed over pairs of N atoms with positions x and radii, r, and where p is self.power and w is the product of self.weight and self.weightFactor Return type:float

Todo

Multiply by the occupancies for simplest approach

Todo

Multiply by sum of occupancies - 1, assuming fully coordinated conformers/rotamers with the simplification that the other conformers do not have atoms that form the same vdW contacts.

VDWRestraint.f_grad(atoms)

Partial derivatives of vdW restraint re: each coordinate & occupancy.

Parameters:atoms (Atoms) – coordinates for current (trial) structure. Returns:Σ:sub:’j|isin|N’ ∂V/∂x:sub:’i’, where V is the van der Waals penalty, for pairs of N atoms with coordinates x. Return type:ndarray[N, 3] QQQ check that not transpose; need to add docc array.

Warning

Does not properly account for open symmetry. Does not work with open symmetry because it ignores certain interactions. If atoms a and b, and their symmetry equivalents a’ and b’, were arranged:

a --- b
b'--- a'

Then motion of a in the y direction would affect the residual of b in a way not taken into account here. Fixing this would greatly increase the computational cost.

Todo

Multiply by the occupancies for simplest approach

Todo

Multiply by sum of occupancies - 1, assuming fully coordinated conformers/rotamers with the simplification that the other conformers do not have atoms that form the same vdW contacts.

Todo

Add partial derivative re: occupancy, and multiply x partials by o. QQQ

VDWRestraint.get_history()

get the history for f and f_grad

VDWRestraint.get_neighbors(atoms, cache=True)

Get the neighbors, using the symmetry module.

Parameters:

• atoms (Atoms) – Unique set of coordinates.
• cache (bool) – try using cached neighbor matrix & radii if change < self.max_shift (always updating the coordinates)

Returns:

symmetry-expanded coordinates, their vdW radii, boolean matrix of which are potentially close

Return type:

Atoms, ndarray, ndarray

Todo

options to exclude same residue, conformer

static VDWRestraint.matrix_bonded(bonds, bondsExcluded=0)

calculate higher degree of bonding connectedness and returns it in matrix form

Parameters:

• bonds (list) – the bonds lists
• bondsExcluded (int) – atom pairs separated by < = bondsExcluded bonds be excluded.

Returns:

matrix with 1 if atoms within bondsExcluded, else 0.

Return type:

numpy.matrix(shape = (n, n), dtype = int)

VDWRestraint.power = 2

Pauli exclusion repels overlapping atoms with force = k.r^11.

This would be too harsh for smooth optimization, so a softer exponent (2) should be used - entirely empirical.

VDWRestraint.radius = {'C': 1.7, 'F': 1.35, 'I': 2.15, 'H': 1.2, 'O': 1.4, 'N': 1.5, 'P': 1.9, 'AS': 2.0, 'BR': 1.95, 'CU': 1.4, 'CL': 1.8, 'S': 1.8, 'SB': 1.85, 'TE': 2.2, 'SE': 2.0}

Default vdW radii.

radiiPath is a hook to enable this to be user over-ridden or extended. However, radiiPath is not yet available as an input argument.

Atoms that are not listed will be given a default radius of 1.5 A.

VDWRestraint.status(atoms)

van der Waals statistics.

Parameters:atoms (Atoms) – coordinates. Returns:object containing number and rms of overlaps, as well as list of worst infractions. Return type:VDWRestraint.Status

VDWRestraint.update(hist, value)

update the history

Parameters:

• hist (str) – ‘f’ or ‘f_grad’ designating the history to extend
• value (str) – text

VDWRestraint.weightFactor = 50

The restraint weight is internally scaled by this, in addition to the user-entered weight. (Avoids user-entered weights being numerically large.)

restraints.ss(atoms, verbose=False)

Identify pairs of Cys sulfur atoms.

Parameters:atoms (Atoms) – coordinates Returns:A_i,j is True when i & j are both SG atoms Return type:ndarray(shape=(natom, natom), dtype=bool)