Author: bugman
Date: Tue Jan 15 20:22:31 2013
New Revision: 18204
URL: http://svn.gna.org/viewcvs/relax?rev=18204&view=rev
Log:
Implemented the pcs.structural_noise user function.
This is used to determine the PCS error due to structural noise via
simulation, and adds the error
via variance addition to the experimental PCS error.
Modified:
trunk/generic_fns/pcs.py
trunk/user_functions/pcs.py
Modified: trunk/generic_fns/pcs.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/generic_fns/pcs.py?rev=18204&r1=18203&r2=18204&view=diff
==============================================================================
--- trunk/generic_fns/pcs.py (original)
+++ trunk/generic_fns/pcs.py Tue Jan 15 20:22:31 2013
@@ -1,6 +1,6 @@
###############################################################################
#
#
-# Copyright (C) 2003-2012 Edward d'Auvergne
#
+# Copyright (C) 2003-2013 Edward d'Auvergne
#
#
#
# This file is part of the program relax (http://www.nmr-relax.com).
#
#
#
@@ -20,13 +20,14 @@
###############################################################################
# Module docstring.
-"""Module for the manipulation of pseudocontact shift data."""
+"""Module for the manipulation of pseudo-contact shift data."""
# Python module imports.
from copy import deepcopy
from math import pi, sqrt
-from numpy import array, float64, ones, zeros
+from numpy import array, float64, ones, std, zeros
from numpy.linalg import norm
+from random import gauss
import sys
from warnings import warn
@@ -34,9 +35,10 @@
from generic_fns import grace, pipes
from generic_fns.align_tensor import get_tensor_index
from generic_fns.mol_res_spin import exists_mol_res_spin_data,
generate_spin_id, return_spin, spin_index_loop, spin_loop
-from maths_fns.pcs import ave_pcs_tensor
+from maths_fns.pcs import ave_pcs_tensor, pcs_tensor
+from maths_fns.vectors import random_unit_vector
from physical_constants import g1H, pcs_constant
-from relax_errors import RelaxError, RelaxAlignError, RelaxNoAlignError,
RelaxNoPdbError, RelaxNoSequenceError, RelaxPCSError
+from relax_errors import RelaxError, RelaxAlignError, RelaxNoAlignError,
RelaxNoPdbError, RelaxNoPCSError, RelaxNoSequenceError, RelaxPCSError
from relax_io import open_write_file, read_spin_data, write_spin_data
from relax_warnings import RelaxWarning, RelaxNoSpinWarning
@@ -666,6 +668,174 @@
cdp.pcs_ids.append(align_id)
+def structural_noise(align_id=None, rmsd=0.2, sim_num=1000, file=None,
dir=None, force=False):
+ """Determine the PCS error due to structural noise via simulation.
+
+ For the simulation the following must already be set up in the current
data pipe:
+
+ - The position of the paramagnetic centre.
+ - The alignment and magnetic susceptibility tensor.
+
+ The protocol for the simulation is as follows:
+
+ - The lanthanide or paramagnetic centre position will be fixed. Its
motion is assumed to be on the femto- to pico- and nanosecond timescales.
Hence the motion is averaged over the evolution of the PCS and can be ignored.
+ - The positions of the nuclear spins will be randomised N times.
For each simulation a random unit vector will be generated. Then a random
distance along the unit vector will be generated by sampling from a Gaussian
distribution centered at zero, the original spin position, with a standard
deviation set to the given RMSD. Both positive and negative displacements
will be used.
+ - The PCS for the randomised position will be back calculated.
+ - The PCS standard deviation will be calculated from the N
randomised PCS values.
+
+ The standard deviation will both be stored in the spin container data
structure in the relax data store as well as being added to the already
present PCS error (using variance addition). This will then be used in any
optimisations involving the PCS.
+
+ If the alignment ID string is not supplied, the procedure will be
applied to the PCS data from all alignments.
+
+
+ @keyword align_id: The alignment tensor ID string.
+ @type align_id: str
+ @keyword rmsd: The atomic position RMSD, in Angstrom, to randomise
the spin positions with for the simulations.
+ @type rmsd: float
+ @keyword sim_num: The number of simulations, N, to perform to
determine the structural noise component of the PCS errors.
+ @type sim_num: int
+ @keyword file: The optional name of the Grace file to plot the
structural errors verses the paramagnetic centre to spin distances.
+ @type file: None or str
+ @keyword dir: The directory name to place the Grace file into.
+ @type dir: None or str
+ @keyword force: A flag which if True will cause any pre-existing
file to be overwritten.
+ @type force: bool
+ """
+
+ # Test if the current pipe exists.
+ pipes.test()
+
+ # Test if sequence data exists.
+ if not exists_mol_res_spin_data():
+ raise RelaxNoSequenceError
+
+ # Test if data corresponding to 'align_id' exists.
+ if not hasattr(cdp, 'pcs_ids'):
+ raise RelaxNoPCSError(align_id)
+ if align_id and align_id not in cdp.pcs_ids:
+ raise RelaxNoPCSError(align_id)
+
+ # Arg check.
+ if align_id and align_id not in cdp.align_ids:
+ raise RelaxError("The alignment ID '%s' is not in the alignment ID
list %s." % (align_id, cdp.align_ids))
+
+ # Convert the align IDs to an array, or take all IDs.
+ if align_id:
+ align_ids = [align_id]
+ else:
+ align_ids = cdp.align_ids
+
+ # Initialise some numpy data structures for use in the simulations.
+ grace_data = []
+ unit_vect = zeros(3, float64)
+ pcs = {}
+ for id in align_ids:
+ grace_data.append([])
+ pcs[id] = zeros(sim_num, float64)
+
+ # Print out.
+ print("Executing %i simulations for each spin system." % sim_num)
+
+ # Loop over the spins.
+ for spin, spin_id in spin_loop(return_id=True):
+ # Deselected spins.
+ if not spin.select:
+ continue
+
+ # Skip spins with no PCS or position.
+ if not hasattr(spin, 'pcs'):
+ continue
+ if not hasattr(spin, 'pos'):
+ continue
+
+ # Print out.
+ print(spin_id)
+
+ # Average the atom position.
+ if type(spin.pos[0]) in [float, float64]:
+ pos = spin.pos
+ else:
+ pos = zeros(3, float64)
+ for i in range(len(spin.pos)):
+ pos += spin.pos[i]
+ pos = pos / len(spin.pos)
+
+ # The original vector length (for the Grace plot).
+ orig_vect = pos - cdp.paramagnetic_centre
+ orig_r = norm(orig_vect)
+
+ # Loop over the N randomisations.
+ for i in range(sim_num):
+ # The random unit vector.
+ random_unit_vector(unit_vect)
+
+ # The random displacement (in Angstrom).
+ disp = gauss(0, rmsd)
+
+ # Move the atom.
+ new_pos = pos + disp*unit_vect
+
+ # The vector and length.
+ vect = new_pos - cdp.paramagnetic_centre
+ r = norm(vect)
+ vect = vect / r
+
+ # Loop over the alignments.
+ for id in align_ids:
+ # No PCS value, so skip.
+ if id not in spin.pcs:
+ continue
+
+ # Calculate the PCS constant.
+ dj = pcs_constant(cdp.temperature[id], cdp.frq[id] * 2.0 *
pi / g1H, r/1e10)
+
+ # Calculate the PCS value (in ppm).
+ pcs[id][i] = pcs_tensor(dj, vect,
cdp.align_tensors[get_tensor_index(id)].A) * 1e6
+
+ # Initialise if necessary.
+ if not hasattr(spin, 'pcs_struct_err'):
+ spin.pcs_struct_err = {}
+
+ # Loop over the alignments.
+ align_index = 0
+ for id in align_ids:
+ # No PCS value, so skip.
+ if id not in spin.pcs:
+ align_index += 1
+ continue
+
+ # The PCS standard deviation.
+ sd = std(pcs[id])
+
+ # Remove the previous error.
+ if id in spin.pcs_struct_err:
+ warn(RelaxWarning("Removing the previous structural error
value from the PCS error of the spin '%s' for the alignment ID '%s'." %
(spin_id, id)))
+ spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 - sd**2)
+
+ # Store the structural error.
+ spin.pcs_struct_err[id] = sd
+
+ # Add it to the PCS error (with variance addition).
+ spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 + sd**2)
+
+ # Store the data for the Grace plot.
+ grace_data[align_index].append([orig_r, sd])
+
+ # Increment the alignment index.
+ align_index += 1
+
+ # The Grace output.
+ if file:
+ # Open the Grace file for writing.
+ file = open_write_file(file, dir, force)
+
+ # The header.
+ grace.write_xy_header(file=file, title="PCS structural noise",
subtitle="%s Angstrom structural noise"%rmsd, sets=len(align_ids),
set_names=align_ids, symbol_sizes=[0.5]*len(align_ids),
linetype=[0]*len(align_ids), data_type=['pcs_bc', 'pcs'],
axis_labels=["Ln\S3+\N to spin distance (Angstrom)", "PCS standard deviation
(ppm)"])
+
+ # The main data.
+ grace.write_xy_data(data=[grace_data], file=file, graph_type='xy')
+
+
def weight(align_id=None, spin_id=None, weight=1.0):
"""Set optimisation weights on the PCS data.
Modified: trunk/user_functions/pcs.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/user_functions/pcs.py?rev=18204&r1=18203&r2=18204&view=diff
==============================================================================
--- trunk/user_functions/pcs.py (original)
+++ trunk/user_functions/pcs.py Tue Jan 15 20:22:31 2013
@@ -1,6 +1,6 @@
###############################################################################
#
#
-# Copyright (C) 2003-2012 Edward d'Auvergne
#
+# Copyright (C) 2003-2013 Edward d'Auvergne
#
#
#
# This file is part of the program relax (http://www.nmr-relax.com).
#
#
#
@@ -391,6 +391,88 @@
uf.wizard_image = WIZARD_IMAGE_PATH + 'align_tensor.png'
+# The pcs.structural_noise user function.
+uf = uf_info.add_uf('pcs.structural_noise')
+uf.title = "Determine the PCS error due to structural noise via simulation."
+uf.title_short = "PCS structural noise simulation."
+uf.display = True
+uf.add_keyarg(
+ name = "align_id",
+ py_type = "str",
+ desc_short = "alignment ID string",
+ desc = "The optional alignment ID string.",
+ wiz_element_type = 'combo',
+ wiz_combo_iter = align_tensor.get_ids,
+ wiz_read_only = True,
+ can_be_none = True
+)
+uf.add_keyarg(
+ name = "rmsd",
+ default = 0.2,
+ py_type = "float",
+ desc_short = "structural RMSD",
+ desc = "The atomic position RMSD, in Angstrom, to randomise the spin
positions with for the simulations."
+)
+uf.add_keyarg(
+ name = "sim_num",
+ default = 1000,
+ min = 3,
+ max = 10000000,
+ py_type = "int",
+ desc_short = "simulation number N",
+ desc = "The number of simulations, N, to perform to determine the
structural noise component of the PCS errors."
+)
+uf.add_keyarg(
+ name = "file",
+ py_type = "str",
+ arg_type = "file sel",
+ desc_short = "Grace file name",
+ desc = "The optional name of the Grace file to plot the structural
errors verses the paramagnetic centre to spin distances.",
+ wiz_filesel_wildcard = "Grace files (*.agr)|*.agr;*.AGR",
+ wiz_filesel_style = FD_SAVE,
+ can_be_none = True
+)
+uf.add_keyarg(
+ name = "dir",
+ py_type = "str",
+ arg_type = "dir",
+ desc_short = "directory name",
+ desc = "The directory name to place the Grace file into.",
+ can_be_none = True
+)
+uf.add_keyarg(
+ name = "force",
+ default = False,
+ py_type = "bool",
+ desc_short = "force flag",
+ desc = "A flag which if True will cause the file to be overwritten."
+)
+# Description.
+uf.desc.append(Desc_container())
+uf.desc[-1].add_paragraph("The analysis of the pseudo-contact shift is
influenced by two significant sources of noise - that of the NMR experiment
and structural noise from the 3D molecular structure used. The closer the
spin to the paramagnetic centre, the greater the influence of structural
noise. This distance dependence is governed by the equation:")
+uf.desc[-1].add_verbatim("""
+ sqrt(3) * abs(delta) * RMSD
+ sigma_dist = --------------------------- ,
+ r
+""")
+uf.desc[-1].add_paragraph("where sigma_dist is the distance component of the
structural noise as a standard deviation, delta is the PCS value, RMSD is the
atomic position root-mean-square deviation, and r is the paramagnetic centre
to spin distance. When close to the paramagnetic centre, this error source
can exceed that of the NMR experiment. The equation for the angular
component of the structural noise is more complicated. The PCS error is
influenced by distance, angle in the alignment frame, and the magnetic
susceptibility tensor.")
+uf.desc[-1].add_paragraph("For the simulation the following must already be
set up in the current data pipe:")
+uf.desc[-1].add_list_element("The position of the paramagnetic centre.")
+uf.desc[-1].add_list_element("The alignment and magnetic susceptibility
tensor.")
+uf.desc[-1].add_paragraph("The protocol for the simulation is as follows:")
+uf.desc[-1].add_list_element("The lanthanide or paramagnetic centre position
will be fixed. Its motion is assumed to be on the femto- to pico- and
nanosecond timescales. Hence the motion is averaged over the evolution of
the PCS and can be ignored.")
+uf.desc[-1].add_list_element("The positions of the nuclear spins will be
randomised N times. For each simulation a random unit vector will be
generated. Then a random distance along the unit vector will be generated by
sampling from a Gaussian distribution centered at zero, the original spin
position, with a standard deviation set to the given RMSD. Both positive and
negative displacements will be used.")
+uf.desc[-1].add_list_element("The PCS for the randomised position will be
back calculated.")
+uf.desc[-1].add_list_element("The PCS standard deviation will be calculated
from the N randomised PCS values.")
+uf.desc[-1].add_paragraph("The standard deviation will both be stored in the
spin container data structure in the relax data store as well as being added
to the already present PCS error (using variance addition). This will then
be used in any optimisations involving the PCS.")
+uf.desc[-1].add_paragraph("If the alignment ID string is not supplied, the
procedure will be applied to the PCS data from all alignments.")
+uf.backend = pcs.structural_noise
+uf.menu_text = "&structural_noise"
+uf.wizard_size = (1000, 700)
+uf.wizard_image = WIZARD_IMAGE_PATH + 'align_tensor.png'
+uf.wizard_apply_button = False
+
+
# The pcs.weight user function.
uf = uf_info.add_uf('pcs.weight')
uf.title = "Set optimisation weights on the PCS data."
r18204 - in /trunk: generic_fns/pcs.py user_functions/pcs.py