Author: bugman
Date: Mon Nov 14 17:00:33 2016
New Revision: 28272
URL: http://svn.gna.org/viewcvs/relax?rev=28272&view=rev
Log:
Bug fix for the pcs.structural_error user function.
The user function now uses a real multivariate normal distribution for
sampling atomic positions.
The previous random unit vector + univariate Gaussian sampling does not
correctly reproduce the
multivariate normal distribution.
Modified:
trunk/pipe_control/pcs.py
trunk/test_suite/system_tests/pcs.py
Modified: trunk/pipe_control/pcs.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/pipe_control/pcs.py?rev=28272&r1=28271&r2=28272&view=diff
==============================================================================
--- trunk/pipe_control/pcs.py (original)
+++ trunk/pipe_control/pcs.py Mon Nov 14 17:00:33 2016
@@ -1,6 +1,6 @@
###############################################################################
#
#
-# Copyright (C) 2003-2015 Edward d'Auvergne
#
+# Copyright (C) 2003-2016 Edward d'Auvergne
#
#
#
# This file is part of the program relax (http://www.nmr-relax.com).
#
#
#
@@ -25,9 +25,9 @@
# Python module imports.
from copy import deepcopy
from math import ceil, floor, pi, sqrt
-from numpy import array, float64, int32, ones, std, zeros
+from numpy import array, eye, float64, int32, ones, std, zeros
from numpy.linalg import norm
-from random import gauss
+from numpy.random import multivariate_normal
import sys
from warnings import warn
@@ -35,7 +35,6 @@
from lib.alignment.pcs import ave_pcs_tensor, pcs_tensor
from lib.check_types import is_float
from lib.errors import RelaxError, RelaxNoAlignError, RelaxNoPdbError,
RelaxNoPCSError, RelaxNoSequenceError
-from lib.geometry.vectors import random_unit_vector
from lib.io import open_write_file, write_data
from lib.periodic_table import periodic_table
from lib.physical_constants import pcs_constant
@@ -1090,7 +1089,7 @@
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 positions of the nuclear spins will be randomised N times
using a multivariate normal distribution.
- The PCS for the randomised position will be back calculated.
- The PCS standard deviation will be calculated from the N
randomised PCS values.
@@ -1130,6 +1129,9 @@
grace_data.append([])
pcs[id] = zeros(sim_num, float64)
+ # Pre-calculate the covariance matrix as the spherical covariance.
+ cov = rmsd**2 * eye(3)
+
# Print out.
print("Executing %i simulations for each spin system." % sim_num)
@@ -1161,19 +1163,13 @@
orig_vect = pos - cdp.paramagnetic_centre
orig_r = norm(orig_vect)
+ # Sample from the multivariate normal distribution.
+ new_pos = multivariate_normal(pos, cov, sim_num)
+
# 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
+ vect = new_pos[i] - cdp.paramagnetic_centre
r = norm(vect)
vect = vect / r
Modified: trunk/test_suite/system_tests/pcs.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/pcs.py?rev=28272&r1=28271&r2=28272&view=diff
==============================================================================
--- trunk/test_suite/system_tests/pcs.py (original)
+++ trunk/test_suite/system_tests/pcs.py Mon Nov 14 17:00:33 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).
#
#
#
@@ -511,16 +511,16 @@
# The simulated data (from 1,000,000 randomisations of 0.2 Angstrom
RMSD).
pcs_struct_err = {
- 'Dy N-dom': 0.014643633242475744,
- 'Er N-dom': 0.0047594540182391868,
- 'Tm N-dom': 0.010454580925459261,
- 'Tb N-dom': 0.01613972832580988
+ 'Dy N-dom': 0.0253,
+ 'Er N-dom': 0.0081,
+ 'Tm N-dom': 0.0181,
+ 'Tb N-dom': 0.0280
}
pcs_err = {
- 'Dy N-dom': 0.1010664929367797,
- 'Er N-dom': 0.10011319794388618,
- 'Tm N-dom': 0.1005450061531003,
- 'Tb N-dom': 0.10129408092495312
+ 'Dy N-dom': 0.1031,
+ 'Er N-dom': 0.1001,
+ 'Tm N-dom': 0.1016,
+ 'Tb N-dom': 0.1038
}
# Alias the single spin.
r28272 - in /trunk: pipe_control/pcs.py test_suite/system_tests/pcs.py