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. _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-commits mailing list relax-commits@gna.org To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-commits