Author: bugman
Date: Thu Jun 2 19:09:39 2016
New Revision: 28222
URL: http://svn.gna.org/viewcvs/relax?rev=28222&view=rev
Log:
Implemented the per-atom RMSD calculation for the structure.rmsd user
function.
Modified:
trunk/lib/structure/statistics.py
trunk/pipe_control/structure/main.py
Modified: trunk/lib/structure/statistics.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/lib/structure/statistics.py?rev=28222&r1=28221&r2=28222&view=diff
==============================================================================
--- trunk/lib/structure/statistics.py (original)
+++ trunk/lib/structure/statistics.py Thu Jun 2 19:09:39 2016
@@ -1,6 +1,6 @@
###############################################################################
#
#
-# Copyright (C) 2011-2015 Edward d'Auvergne
#
+# Copyright (C) 2011-2016 Edward d'Auvergne
#
#
#
# This file is part of the program relax (http://www.nmr-relax.com).
#
#
#
@@ -115,3 +115,42 @@
# Average.
mean[i] = mean[i] / weights.sum()
+
+
+def per_atom_rmsd(coord, verbosity=0):
+ """Determine the per-atom RMSDs 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
+ @return: The list of RMSD values for each atom.
+ @rtype: rank-1 numpy float64 array
+ """
+
+ # Init.
+ M = len(coord)
+ N = len(coord[0])
+ model_rmsd = zeros(M, float64)
+ mean_str = zeros((N, 3), float64)
+ rmsd = zeros(N, float64)
+
+ # Calculate the mean structure.
+ calc_mean_structure(coord, mean_str)
+
+ # Loop over the atoms.
+ for j in range(N):
+ # Loop over the models.
+ for i in range(M):
+ # The vector connecting the mean to model atom.
+ vect = mean_str[j] - coord[i][j]
+
+ # The atomic RMSD.
+ rmsd[j] += norm(vect)**2
+
+ # Normalise, and sqrt.
+ rmsd[j] = sqrt(rmsd[j] / M)
+
+ # Return the RMSDs.
+ return rmsd
Modified: trunk/pipe_control/structure/main.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/pipe_control/structure/main.py?rev=28222&r1=28221&r2=28222&view=diff
==============================================================================
--- trunk/pipe_control/structure/main.py (original)
+++ trunk/pipe_control/structure/main.py Thu Jun 2 19:09:39 2016
@@ -44,7 +44,7 @@
from lib.structure.internal.object import Internal
from lib.structure.pca import pca_analysis
from lib.structure.represent.diffusion_tensor import diffusion_tensor
-from lib.structure.statistics import atomic_rmsd
+from lib.structure.statistics import atomic_rmsd, per_atom_rmsd
from lib.structure.superimpose import fit_to_first, fit_to_mean
from lib.warnings import RelaxWarning, RelaxNoPDBFileWarning,
RelaxZeroVectorWarning
from pipe_control import molmol, pipes
@@ -1308,7 +1308,7 @@
cdp.structure.load_xyz(file_path, read_mol=read_mol,
set_mol_name=set_mol_name, read_model=read_model,
set_model_num=set_model_num, verbosity=verbosity)
-def rmsd(pipes=None, models=None, molecules=None, atom_id=None):
+def rmsd(pipes=None, models=None, molecules=None, atom_id=None,
atomic=False):
"""Calculate the RMSD between the loaded models.
The RMSD value will be placed into the current data pipe
cdp.structure.rmsd data structure.
@@ -1322,6 +1322,8 @@
@type molecules: None or list of lists of str
@keyword atom_id: The atom identification string of the coordinates of
interest. This matches the spin ID string format.
@type atom_id: str or None
+ @keyword atomic: A flag which if True will allow for per-atom RMSDs
to be additionally calculated.
+ @type atomic: bool
@return: The RMSD value.
@rtype: float
"""
@@ -1332,10 +1334,34 @@
# Assemble the structural coordinates.
coord, ids, mol_names, res_names, res_nums, atom_names, elements =
assemble_structural_coordinates(pipes=pipes, models=models,
molecules=molecules, atom_id=atom_id)
+ # Per-atom RMSDs.
+ if atomic:
+ # Printout.
+ print("\nCalculating the atomic-level RMSDs.")
+
+ # Calculate the per-atom RMSDs.
+ rmsd = per_atom_rmsd(coord, verbosity=0)
+
+ # Loop over the atoms.
+ for i in range(len(res_nums)):
+ # The spin identification string.
+ id = generate_spin_id_unique(mol_name=mol_names[i],
res_num=res_nums[i], res_name=res_names[i], spin_num=i,
spin_name=atom_names[i])
+
+ # Get the spin container.
+ spin_cont = return_spin(id)
+
+ # Skip the spin if it doesn't exist.
+ if spin_cont == None:
+ continue
+
+ # Store the value.
+ spin_cont.pos_rmsd = rmsd[i]
+
# Calculate the RMSD.
+ print("\nCalculating the global RMSD.")
cdp.structure.rmsd = atomic_rmsd(coord, verbosity=1)
- # Return the RMSD.
+ # Return the global RMSD.
return cdp.structure.rmsd
r28222 - in /trunk: lib/structure/statistics.py pipe_control/structure/main.py