Author: bugman
Date: Thu Nov 24 14:04:42 2011
New Revision: 15012
URL: http://svn.gna.org/viewcvs/relax?rev=15012&view=rev
Log:
Created a new algorithm for finding the pivot of motion between different
structural models.
This is available through the structure.find_pivot user function.
Modified:
1.3/generic_fns/structure/main.py
1.3/generic_fns/structure/superimpose.py
1.3/prompt/structure.py
Modified: 1.3/generic_fns/structure/main.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/structure/main.py?rev=15012&r1=15011&r2=15012&view=diff
==============================================================================
--- 1.3/generic_fns/structure/main.py (original)
+++ 1.3/generic_fns/structure/main.py Thu Nov 24 14:04:42 2011
@@ -22,6 +22,7 @@
# Python module imports.
from math import sqrt
+from minfx.generic import generic_minimise
from numpy import array, dot, float64, ndarray, zeros
from numpy.linalg import norm
from os import F_OK, access
@@ -37,7 +38,7 @@
from generic_fns.structure.api_base import Displacements
from generic_fns.structure.internal import Internal
from generic_fns.structure.scientific import Scientific_data
-from generic_fns.structure.superimpose import fit_to_first, fit_to_mean
+from generic_fns.structure.superimpose import fit_to_first, fit_to_mean,
Pivot_finder
from relax_errors import RelaxError, RelaxFileError, RelaxNoPdbError,
RelaxNoSequenceError
from relax_io import get_file_path, open_write_file, write_spin_data
from relax_warnings import RelaxWarning, RelaxNoPDBFileWarning,
RelaxZeroVectorWarning
@@ -168,6 +169,57 @@
# Send to the base container for the calculations.
cdp.structure.displacements._calculate(model_from=model_from[i],
model_to=model_to[j], coord_from=array(coord_from), coord_to=array(coord_to),
centroid=centroid)
+
+
+def find_pivot(models=None, atom_id=None, init_pos=None):
+ """Superimpose a set of structural models.
+
+ @keyword models: The list of models to use. If set to None, then all
models will be used.
+ @type models: list of int or None
+ @keyword atom_id: The molecule, residue, and atom identifier string.
This matches the spin ID string format.
+ @type atom_id: str or None
+ @keyword init_pos: The starting pivot position for the pivot point
optimisation.
+ @type init_pos: list of float or numpy rank-1, 3D array
+ """
+
+ # Initialised the starting position if needed.
+ if init_pos == None:
+ init_pos = zeros(3, float64)
+ init_pos = array(init_pos)
+
+ # Validate the models.
+ cdp.structure.validate_models()
+
+ # Create a list of all models.
+ if models == None:
+ models = []
+ for model in cdp.structure.model_loop():
+ models.append(model.num)
+
+ # Assemble the atomic coordinates of all models.
+ coord = []
+ for model in models:
+ coord.append([])
+ for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model,
pos_flag=True):
+ coord[-1].append(pos[0])
+ coord[-1] = array(coord[-1])
+ coord = array(coord)
+
+ # The target function.
+ finder = Pivot_finder(models, coord)
+ results = generic_minimise(func=finder.func, x0=init_pos,
min_algor='simplex', print_flag=1)
+
+ # No result.
+ if results == None:
+ return
+
+ # Store the data.
+ cdp.structure.pivot = results
+
+ # Print out.
+ print("Motional pivot found at: %s" % results)
+
+
def get_pos(spin_id=None, str_id=None, ave_pos=False):
Modified: 1.3/generic_fns/structure/superimpose.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/structure/superimpose.py?rev=15012&r1=15011&r2=15012&view=diff
==============================================================================
--- 1.3/generic_fns/structure/superimpose.py (original)
+++ 1.3/generic_fns/structure/superimpose.py Thu Nov 24 14:04:42 2011
@@ -26,11 +26,90 @@
# Python module imports.
from copy import deepcopy
from math import pi
-from numpy import diag, dot, eye, float64, outer, sign, transpose, zeros
+from numpy import diag, dot, eye, float64, outer, sign, sqrt, transpose,
zeros
from numpy.linalg import det, norm, svd
# relax module import.
from maths_fns.rotation_matrix import R_to_axis_angle
+
+
+class Pivot_finder:
+ """Class for finding the optimal pivot point for motions between the
given models."""
+
+ def __init__(self, models, coord):
+ """Set up the class for pivot point optimisation.
+
+ @keyword models: The list of models to use. If set to None, then
all models will be used.
+ @type models: list of int or None
+ @keyword coord: The array of molecular coordinates. The first
dimension corresponds to the model, the second the atom, the third the
coordinate.
+ @type coord: rank-3 numpy array
+ """
+
+ # Store the args.
+ self.models = models
+ self.coord = coord
+
+ # Store a copy of the coordinates for restoration.
+ self.orig_coord = deepcopy(coord)
+
+
+ def func(self, params):
+ """Target function for the optimisation of the motional pivot point.
+
+ @param params: The parameter vector from the optimisation algorithm.
+ @type params: list
+ @return: The target function value defined as the combined
RMSD value.
+ @rtype: float
+ """
+
+ # The fit to mean algorithm.
+ T, R, pivot = fit_to_mean(models=self.models, coord=self.coord,
centroid=params, verbosity=0)
+
+ # The RMSD.
+ val = atomic_rmsd(self.coord)
+
+ # Restore the coordinates.
+ self.coord = deepcopy(self.orig_coord)
+
+ # Return the RMSD.
+ return val
+
+
+
+def atomic_rmsd(coord):
+ """Determine the RMSD for the given atomic coordinates.
+
+ This is the per atom RMSD to the mean structure.
+
+
+ @keyword coord: The array of molecular coordinates. The first
dimension corresponds to the model, the second the atom, the third the
coordinate.
+ @type coord: rank-3 numpy array
+ """
+
+ # Init.
+ M = len(coord)
+ N = len(coord[0])
+ model_rmsd = zeros(M, float64)
+ mean = zeros((N, 3), float64)
+
+ # Calculate the mean structure.
+ calc_mean_structure(coord, mean)
+
+ # Loop over the models.
+ for i in range(M):
+ # Loop over the atoms.
+ for j in range(N):
+ # The vector connecting the mean to model atom.
+ vect = mean[j] - coord[i][j]
+
+ # The atomic RMSD.
+ model_rmsd[i] += norm(vect)**2
+
+ # Normalise, and sqrt.
+ model_rmsd[i] = sqrt(model_rmsd[i] / N)
+
+ # Return the average RMSD.
+ return sum(model_rmsd) / M
def calc_mean_structure(coord=None, mean=None):
@@ -111,7 +190,7 @@
return T_list, R_list, pivot_list
-def fit_to_mean(models=None, coord=None, centroid=None):
+def fit_to_mean(models=None, coord=None, centroid=None, verbosity=1):
"""Superimpose a set of structural models using the fit to first
algorithm.
@keyword models: The list of models to superimpose.
@@ -120,12 +199,15 @@
@type coord: list of numpy rank-2, Nx3 arrays
@keyword centroid: An alternative position of the centroid to allow for
different superpositions, for example of pivot point motions.
@type centroid: list of float or numpy rank-1, 3D array
+ @keyword verbosity: The amount of information to print out. If 0,
nothing will be printed.
+ @type verbosity: int
@return: The lists of translation vectors, rotation matrices,
and rotation pivots.
@rtype: list of numpy rank-1 3D arrays, list of numpy rank-2
3D arrays, list of numpy rank-1 3D arrays
"""
# Print out.
- print("\nSuperimposition of structural models %s using the 'fit to mean'
algorithm." % models)
+ if verbosity:
+ print("\nSuperimposition of structural models %s using the 'fit to
mean' algorithm." % models)
# Duplicate the coordinates.
orig_coord = deepcopy(coord)
@@ -144,8 +226,9 @@
iter = 0
while not converged:
# Print out.
- print("\nIteration %i of the algorithm." % iter)
- print("%-10s%-25s%-25s" % ("Model", "Translation (Angstrom)",
"Rotation (deg)"))
+ if verbosity:
+ print("\nIteration %i of the algorithm." % iter)
+ print("%-10s%-25s%-25s" % ("Model", "Translation (Angstrom)",
"Rotation (deg)"))
# Calculate the mean structure.
calc_mean_structure(coord, mean)
@@ -157,10 +240,10 @@
trans_vect, trans_dist, R, axis, angle, pivot =
kabsch(name_from='model %s'%models[0], name_to='mean', coord_from=coord[i],
coord_to=mean, centroid=centroid, verbosity=0)
# Table print out.
- print("%-10i%25.3g%25.3g" % (i, trans_dist, (angle / 2.0 / pi *
360.0)))
+ if verbosity:
+ print("%-10i%25.3g%25.3g" % (i, trans_dist, (angle / 2.0 /
pi * 360.0)))
# Shift the coordinates.
- print coord[i][0]
for j in range(N):
# Translate.
coord[i][j] += trans_vect
@@ -173,7 +256,6 @@
# The new position.
coord[i][j] += pivot
- print coord[i][0]
# Convergence test.
if trans_dist > 1e-10 or angle > 1e-10:
@@ -185,7 +267,7 @@
# Perform the fit once from the original coordinates to obtain the full
transforms.
for i in range(len(models)):
# Calculate the displacements (Kabsch algorithm).
- trans_vect, trans_dist, R, axis, angle, pivot =
kabsch(name_from='model %s'%models[i], name_to='the mean structure',
coord_from=orig_coord[i], coord_to=mean, centroid=centroid)
+ trans_vect, trans_dist, R, axis, angle, pivot =
kabsch(name_from='model %s'%models[i], name_to='the mean structure',
coord_from=orig_coord[i], coord_to=mean, centroid=centroid, verbosity=0)
# Store the transforms.
T_list.append(trans_vect)
Modified: 1.3/prompt/structure.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/prompt/structure.py?rev=15012&r1=15011&r2=15012&view=diff
==============================================================================
--- 1.3/prompt/structure.py (original)
+++ 1.3/prompt/structure.py Thu Nov 24 14:04:42 2011
@@ -332,6 +332,41 @@
_build_doc(displacement)
+ def find_pivot(self, models=None, atom_id=None, init_pos=None):
+ # Function intro text.
+ if self._exec_info.intro:
+ text = self._exec_info.ps3 + "structure.find_pivot("
+ text = text + "models=" + repr(models)
+ text = text + ", atom_id=" + repr(atom_id)
+ text = text + ", init_pos=" + repr(init_pos) + ")"
+ print(text)
+
+ # The argument checks.
+ arg_check.is_int_list(models, 'model list', can_be_none=True)
+ arg_check.is_str(atom_id, 'atom identification string',
can_be_none=True)
+ arg_check.is_float_array(init_pos, 'initial pivot position',
can_be_none=True)
+
+ # Execute the functional code.
+ generic_fns.structure.main.find_pivot(models=models,
atom_id=atom_id, init_pos=init_pos)
+
+ # The function doc info.
+ find_pivot._doc_title = "Find the pivot point of the motion of a set of
structures."
+ find_pivot._doc_title_short = "Pivot search."
+ find_pivot._doc_args = [
+ ["models", "The list of models to use."],
+ ["atom_id", "The atom identification string."],
+ ["init_pos", "The initial position of the pivot."]
+ ]
+ find_pivot._doc_desc = """
+ This is used to find pivot point of motion between a set of
structural models. If the list of models is not supplied, then all models
will be used.
+
+ The atom ID, which uses the same notation as the spin ID strings,
can be used to restrict the search to certain molecules, residues, or atoms.
For example to only use backbone heavy atoms in a protein, use the atom ID of
'@N,C,CA,O', assuming those are the names of the atoms from the structural
file.
+
+ By supplying the position of the centroid, an alternative position
than the standard rigid body centre is used as the focal point of the
superimposition. The allows, for example, the superimposition about a pivot
point.
+ """
+ _build_doc(find_pivot)
+
+
def load_spins(self, spin_id=None, ave_pos=True):
# Function intro text.
if self._exec_info.intro:
r15012 - in /1.3: generic_fns/structure/main.py generic_fns/structure/superimpose.py prompt/structure.py