Author: bugman Date: Mon Feb 3 13:35:23 2014 New Revision: 22110 URL: http://svn.gna.org/viewcvs/relax?rev=22110&view=rev Log: Initial implementation of the lib.frame_order.double_rotor module. This module implements the functions needed to solve the frame order analysis for the RDC (via the frame order matrix) and PCS (numerically). The interfaces have been updated for the double rotor but most of the code still implements the basic rotor model from which it derives. Added: branches/double_rotor/lib/frame_order/double_rotor.py - copied, changed from r22105, branches/double_rotor/lib/frame_order/rotor.py Modified: branches/double_rotor/lib/frame_order/__init__.py branches/double_rotor/target_functions/frame_order.py Modified: branches/double_rotor/lib/frame_order/__init__.py URL: http://svn.gna.org/viewcvs/relax/branches/double_rotor/lib/frame_order/__init__.py?rev=22110&r1=22109&r2=22110&view=diff ============================================================================== --- branches/double_rotor/lib/frame_order/__init__.py (original) +++ branches/double_rotor/lib/frame_order/__init__.py Mon Feb 3 13:35:23 2014 @@ -1,6 +1,6 @@ ############################################################################### # # -# Copyright (C) 2013 Edward d'Auvergne # +# Copyright (C) 2013-2014 Edward d'Auvergne # # # # This file is part of the program relax (http://www.nmr-relax.com). # # # @@ -23,6 +23,7 @@ """The relax-lib NMR package - a library of functions for the frame order theories.""" __all__ = [ + 'double_rotor', 'format', 'free_rotor', 'iso_cone_free_rotor', Copied: branches/double_rotor/lib/frame_order/double_rotor.py (from r22105, branches/double_rotor/lib/frame_order/rotor.py) URL: http://svn.gna.org/viewcvs/relax/branches/double_rotor/lib/frame_order/double_rotor.py?p2=branches/double_rotor/lib/frame_order/double_rotor.py&p1=branches/double_rotor/lib/frame_order/rotor.py&r1=22105&r2=22110&rev=22110&view=diff ============================================================================== --- branches/double_rotor/lib/frame_order/rotor.py (original) +++ branches/double_rotor/lib/frame_order/double_rotor.py Mon Feb 3 13:35:23 2014 @@ -1,6 +1,6 @@ ############################################################################### # # -# Copyright (C) 2009-2013 Edward d'Auvergne # +# Copyright (C) 2009-2014 Edward d'Auvergne # # # # This file is part of the program relax (http://www.nmr-relax.com). # # # @@ -20,7 +20,7 @@ ############################################################################### # Module docstring. -"""Module for the handling of Frame Order.""" +"""Module for the double rotor frame order model.""" # Dependency check module. import dep_check @@ -36,8 +36,8 @@ from lib.frame_order.matrix_ops import rotate_daeg -def compile_2nd_matrix_rotor(matrix, Rx2_eigen, smax): - """Generate the rotated 2nd degree Frame Order matrix for the rotor model. +def compile_2nd_matrix_double_rotor(matrix, Rx2_eigen, smax, smaxb): + """Generate the rotated 2nd degree Frame Order matrix for the double rotor model. The cone axis is assumed to be parallel to the z-axis in the eigenframe. @@ -46,8 +46,10 @@ @type matrix: numpy 9D, rank-2 array @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself. @type Rx2_eigen: numpy 9D, rank-2 array - @param smax: The maximum torsion angle. + @param smax: The maximum torsion angle for the first rotor. @type smax: float + @param smaxb: The maximum torsion angle for the second rotor. + @type smaxb: float """ # Zeros. @@ -80,53 +82,15 @@ return rotate_daeg(matrix, Rx2_eigen) -def pcs_numeric_int_rotor(sigma_max=None, c=None, r_pivot_atom=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None): - """Determine the averaged PCS value via numerical integration. - - @keyword sigma_max: The maximum rotor angle. - @type sigma_max: float - @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). - @type c: float - @keyword r_pivot_atom: The pivot point to atom vector. - @type r_pivot_atom: numpy rank-1, 3D array - @keyword r_ln_pivot: The lanthanide position to pivot point vector. - @type r_ln_pivot: numpy rank-1, 3D array - @keyword A: The full alignment tensor of the non-moving domain. - @type A: numpy rank-2, 3D array - @keyword R_eigen: The eigenframe rotation matrix. - @type R_eigen: numpy rank-2, 3D array - @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). - @type RT_eigen: numpy rank-2, 3D array - @keyword Ri_prime: The empty rotation matrix for the in-frame rotor motion, used to calculate the PCS for each state i in the numerical integration. - @type Ri_prime: numpy rank-2, 3D array - @return: The averaged PCS value. - @rtype: float - """ - - # Preset the rotation matrix elements for state i. - Ri_prime[0, 2] = 0.0 - Ri_prime[1, 2] = 0.0 - Ri_prime[2, 0] = 0.0 - Ri_prime[2, 1] = 0.0 - Ri_prime[2, 2] = 1.0 - - # Perform numerical integration. - result = quad(pcs_pivot_motion_rotor, -sigma_max, sigma_max, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime)) - - # The surface area normalisation factor. - SA = 2.0 * sigma_max - - # Return the value. - return c * result[0] / SA - - -def pcs_numeric_int_rotor_qrint(points=None, sigma_max=None, c=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False): - """Determine the averaged PCS value via numerical integration. +def pcs_numeric_int_double_rotor(points=None, sigma_max=None, sigma_max_2=None, c=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False): + """The averaged PCS value via numerical integration for the double rotor frame order model. @keyword points: The Sobol points in the torsion-tilt angle space. @type points: numpy rank-2, 3D array - @keyword sigma_max: The maximum rotor angle. + @keyword sigma_max: The maximum opening angle for the first rotor. @type sigma_max: float + @keyword sigma_max_2: The maximum opening angle for the second rotor. + @type sigma_max_2: float @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). @type c: numpy rank-1 array @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. @@ -165,14 +129,14 @@ num = 0 for i in range(len(points)): # Unpack the point. - sigma = points[i] + sigma, sigma2 = points[i] # Outside of the distribution, so skip the point. if sigma > sigma_max or sigma < -sigma_max: continue # Calculate the PCSs for this state. - pcs_pivot_motion_rotor_qrint(sigma_i=sigma, full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom, r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A, R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs) + pcs_pivot_motion_double_rotor(sigma_i=sigma, full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom, r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A, R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs) # Increment the number of points. num += 1 @@ -190,56 +154,8 @@ print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6)) -def pcs_pivot_motion_rotor(sigma_i, r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime): - """Calculate the PCS value after a pivoted motion for the rotor model. - - @param sigma_i: The rotor angle for state i. - @type sigma_i: float - @param r_pivot_atom: The pivot point to atom vector. - @type r_pivot_atom: numpy rank-1, 3D array - @param r_ln_pivot: The lanthanide position to pivot point vector. - @type r_ln_pivot: numpy rank-1, 3D array - @param A: The full alignment tensor of the non-moving domain. - @type A: numpy rank-2, 3D array - @param R_eigen: The eigenframe rotation matrix. - @type R_eigen: numpy rank-2, 3D array - @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). - @type RT_eigen: numpy rank-2, 3D array - @param Ri_prime: The empty rotation matrix for the in-frame rotor motion for state i. - @type Ri_prime: numpy rank-2, 3D array - @return: The PCS value for the changed position. - @rtype: float - """ - - # The rotation matrix. - c_sigma = cos(sigma_i) - s_sigma = sin(sigma_i) - Ri_prime[0, 0] = c_sigma - Ri_prime[0, 1] = -s_sigma - Ri_prime[1, 0] = s_sigma - Ri_prime[1, 1] = c_sigma - - # The rotation. - R_i = dot(R_eigen, dot(Ri_prime, RT_eigen)) - - # Calculate the new vector. - vect = dot(R_i, r_pivot_atom) + r_ln_pivot - - # The vector length. - length = norm(vect) - - # The projection. - proj = dot(vect, dot(A, vect)) - - # The PCS. - pcs = proj / length**5 - - # Return the PCS value (without the PCS constant). - return pcs - - -def pcs_pivot_motion_rotor_qrint(sigma_i=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False): - """Calculate the PCS value after a pivoted motion for the rotor model. +def pcs_pivot_motion_double_rotor(sigma_i=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False): + """Calculate the PCS value after a pivoted motion for the double rotor model. @keyword sigma_i: The rotor angle for state i. @type sigma_i: float Modified: branches/double_rotor/target_functions/frame_order.py URL: http://svn.gna.org/viewcvs/relax/branches/double_rotor/target_functions/frame_order.py?rev=22110&r1=22109&r2=22110&view=diff ============================================================================== --- branches/double_rotor/target_functions/frame_order.py (original) +++ branches/double_rotor/target_functions/frame_order.py Mon Feb 3 13:35:23 2014 @@ -35,6 +35,7 @@ from lib.alignment.rdc import rdc_tensor from lib.errors import RelaxError from lib.float import isNaN +from lib.frame_order.double_rotor import compile_2nd_matrix_double_rotor, pcs_numeric_int_double_rotor from lib.frame_order.free_rotor import compile_2nd_matrix_free_rotor from lib.frame_order.iso_cone import compile_2nd_matrix_iso_cone, pcs_numeric_int_iso_cone, pcs_numeric_int_iso_cone_qrint from lib.frame_order.iso_cone_free_rotor import compile_2nd_matrix_iso_cone_free_rotor