Author: bugman
Date: Tue Mar 25 17:36:10 2014
New Revision: 22529
URL: http://svn.gna.org/viewcvs/relax?rev=22529&view=rev
Log:
Shifted all of the code for calculating the frame order rotor axis into
lib.frame_order.rotor_axis.
The new frame_order.rotor_axis module consists of three function for creating
a unit vector or the
rotor axis using either the axis alpha angle, the two spherical angles or the
three Euler angles.
Added:
branches/double_rotor/lib/frame_order/rotor_axis.py
Modified:
branches/double_rotor/lib/frame_order/__init__.py
branches/double_rotor/specific_analyses/frame_order/user_functions.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=22529&r1=22528&r2=22529&view=diff
==============================================================================
--- branches/double_rotor/lib/frame_order/__init__.py (original)
+++ branches/double_rotor/lib/frame_order/__init__.py Tue Mar 25 17:36:10
2014
@@ -33,5 +33,6 @@
'pseudo_ellipse_free_rotor',
'pseudo_ellipse',
'pseudo_ellipse_torsionless',
- 'rotor'
+ 'rotor',
+ 'rotor_axis'
]
Added: branches/double_rotor/lib/frame_order/rotor_axis.py
URL:
http://svn.gna.org/viewcvs/relax/branches/double_rotor/lib/frame_order/rotor_axis.py?rev=22529&view=auto
==============================================================================
--- branches/double_rotor/lib/frame_order/rotor_axis.py (added)
+++ branches/double_rotor/lib/frame_order/rotor_axis.py Tue Mar 25 17:36:10
2014
@@ -0,0 +1,108 @@
+###############################################################################
+#
#
+# Copyright (C) 2014 Edward d'Auvergne
#
+#
#
+# This file is part of the program relax (http://www.nmr-relax.com).
#
+#
#
+# This program is free software: you can redistribute it and/or modify
#
+# it under the terms of the GNU General Public License as published by
#
+# the Free Software Foundation, either version 3 of the License, or
#
+# (at your option) any later version.
#
+#
#
+# This program is distributed in the hope that it will be useful,
#
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
#
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#
+# GNU General Public License for more details.
#
+#
#
+# You should have received a copy of the GNU General Public License
#
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
+#
#
+###############################################################################
+
+# Module docstring.
+"""Functions for creating or calculating the rotor axis for the frame order
models."""
+
+# Python module imports.
+from numpy import array, cross, dot, float64, zeros
+from numpy.linalg import norm
+
+# relax module imports.
+from lib.geometry.coord_transform import spherical_to_cartesian
+from lib.geometry.rotations import axis_angle_to_R, euler_to_R_zyz
+
+# Module variables.
+R = zeros((3, 3), float64) # A rotation matrix.
+Z_AXIS = array([0, 0, 1], float64)
+
+
+def create_rotor_axis_alpha(alpha=None, pivot=None, point=None):
+ """Create the rotor axis from the axis alpha angle.
+
+ @keyword alpha: The axis alpha angle, defined as the angle between a
vector perpendicular to the pivot-CoM vector in the xy-plane and the rotor
axis.
+ @type alpha: float
+ @keyword pivot: The pivot point on the rotation axis.
+ @type pivot: numpy rank-1 3D array
+ @keyword point: The reference point in space.
+ @type point: numpy rank-1 3D array
+ @return: The rotor axis as a unit vector.
+ @rtype: numpy rank-1 3D float64 array
+ """
+
+ # The point-pivot unit vector - the norm of the system (the pivot is
defined as the point on the axis closest to the point).
+ n = point - pivot
+ n = n / norm(n)
+
+ # The vector perpendicular to the CoM-pivot vector and in the xy plane.
+ mu_xy = cross(n, Z_AXIS)
+ mu_xy = mu_xy / norm(mu_xy)
+
+ # Rotate the vector about the CoM-pivot axis by the angle alpha.
+ axis_angle_to_R(n, alpha, R)
+ axis = dot(R, mu_xy)
+
+ # Return the new axis.
+ return axis
+
+
+def create_rotor_axis_euler(alpha=None, beta=None, gamma=None):
+ """Create the rotor axis from the Euler angles.
+
+ @keyword alpha: The alpha Euler angle in the zyz notation.
+ @type alpha: float
+ @keyword beta: The beta Euler angle in the zyz notation.
+ @type beta: float
+ @keyword gamma: The gamma Euler angle in the zyz notation.
+ @type gamma: float
+ @return: The rotor axis as a unit vector.
+ @rtype: numpy rank-1 3D float64 array
+ """
+
+ # Initialise the 3D frame.
+ frame = zeros((3, 3), float64)
+
+ # Euler angle to rotation matrix conversion.
+ euler_to_R_zyz(alpha, beta, gamma, frame)
+
+ # Return the z-axis component.
+ return frame[:, 2]
+
+
+def create_rotor_axis_spherical(theta=None, phi=None):
+ """Create the rotor axis from the spherical coordinates.
+
+ @keyword theta: The polar spherical angle.
+ @type theta: float
+ @keyword phi: The azimuthal spherical angle.
+ @type phi: float
+ @return: The rotor axis as a unit vector.
+ @rtype: numpy rank-1 3D float64 array
+ """
+
+ # Initialise the axis.
+ axis = zeros(3, float64)
+
+ # Parameter conversion.
+ spherical_to_cartesian([1.0, theta, phi], axis)
+
+ # Return the new axis.
+ return axis
Modified:
branches/double_rotor/specific_analyses/frame_order/user_functions.py
URL:
http://svn.gna.org/viewcvs/relax/branches/double_rotor/specific_analyses/frame_order/user_functions.py?rev=22529&r1=22528&r2=22529&view=diff
==============================================================================
--- branches/double_rotor/specific_analyses/frame_order/user_functions.py
(original)
+++ branches/double_rotor/specific_analyses/frame_order/user_functions.py
Tue Mar 25 17:36:10 2014
@@ -31,7 +31,7 @@
# relax module imports.
from lib.arg_check import is_float_array
from lib.errors import RelaxError
-from lib.geometry.coord_transform import spherical_to_cartesian
+from lib.frame_order.rotor_axis import create_rotor_axis_alpha,
create_rotor_axis_euler, create_rotor_axis_spherical
from lib.geometry.rotations import euler_to_R_zyz, two_vect_to_R
from lib.io import open_write_file
from lib.order import order_parameters
@@ -181,17 +181,16 @@
else:
rotor_angle = cdp.cone_sigma_max
- # Generate the rotor axis.
- if cdp.model in ['rotor', 'free rotor', 'iso cone', 'iso cone, free
rotor']:
- axis = zeros(3, float64)
- spherical_to_cartesian([1.0, cdp.axis_theta, cdp.axis_phi], axis)
- else:
- axes = zeros((3, 3), float64)
- euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma,
axes)
- axis = axes[:, 2]
-
# Get the CoM of the entire molecule to use as the centre of the
rotor.
com = pipe_centre_of_mass(verbosity=0)
+
+ # Generate the rotor axis.
+ if cdp.model in ['rotor']:
+ axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha,
pivot=cdp.pivot, point=com)
+ elif cdp.model in ['free rotor', 'iso cone', 'iso cone, free rotor']:
+ axis = create_rotor_axis_spherical(theta=cdp.axis_theta,
phi=cdp.axis_phi)
+ else:
+ axis = create_rotor_axis_euler(alpha=cdp.eigen_alpha,
beta=cdp.eigen_beta, gamma=cdp.eigen_gamma)
# Add the rotor object to the structure as a new molecule.
rotor_pdb(structure=structure, rotor_angle=rotor_angle, axis=axis,
axis_pt=cdp.pivot, centre=com, span=2e-9, blade_length=5e-10, staggered=False)
@@ -228,8 +227,7 @@
print("\nGenerating the z-axis system.")
# The axis.
- axis = zeros(3, float64)
- spherical_to_cartesian([1.0, getattr(cdp, 'axis_theta'),
getattr(cdp, 'axis_phi')], axis)
+ axis = create_rotor_axis_spherical(theta=cdp.axis_theta,
phi=cdp.axis_phi)
print(("Central axis: %s." % axis))
# Rotations and inversions.
@@ -245,7 +243,7 @@
# Fill the structure.
for i in range(cdp.sim_number):
- spherical_to_cartesian([1.0, getattr(cdp,
'axis_theta_sim')[i], getattr(cdp, 'axis_phi_sim')[i]], axis_sim[i])
+ axis_sim[i] =
create_rotor_axis_spherical(theta=cdp.axis_theta_sim[i],
phi=cdp.axis_phi_sim[i])
# Inversion.
axis_sim_pos = axis_sim
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=22529&r1=22528&r2=22529&view=diff
==============================================================================
--- branches/double_rotor/target_functions/frame_order.py (original)
+++ branches/double_rotor/target_functions/frame_order.py Tue Mar 25
17:36:10 2014
@@ -45,6 +45,7 @@
from lib.frame_order.pseudo_ellipse_free_rotor import
compile_2nd_matrix_pseudo_ellipse_free_rotor
from lib.frame_order.pseudo_ellipse_torsionless import
compile_2nd_matrix_pseudo_ellipse_torsionless,
pcs_numeric_int_pseudo_ellipse_torsionless,
pcs_numeric_int_pseudo_ellipse_torsionless_qrint
from lib.frame_order.rotor import compile_2nd_matrix_rotor,
pcs_numeric_int_rotor, pcs_numeric_int_rotor_qrint
+from lib.frame_order.rotor_axis import create_rotor_axis_alpha
from lib.geometry.coord_transform import spherical_to_cartesian
from lib.geometry.rotations import axis_angle_to_R, euler_to_R_zyz,
two_vect_to_R
from lib.geometry.vectors import vector_angle
@@ -1896,17 +1897,8 @@
else:
ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha,
sigma_max = params
- # The CoM-pivot unit vector (the pivot is defined as the point on
the axis closest to the CoM).
- piv_com = self.com - array(self._param_pivot, float64)
- piv_com = piv_com / norm(piv_com)
-
- # The vector perpendicular to the CoM-pivot vector and in the xy
plane.
- mu_xy = cross(piv_com, self.z_axis)
- mu_xy = mu_xy / norm(mu_xy)
-
- # Rotate the vector about the CoM-pivot axis by the angle alpha.
- axis_angle_to_R(piv_com, axis_alpha, self.R)
- self.cone_axis = dot(self.R, mu_xy)
+ # Generate the rotor axis.
+ self.cone_axis = create_rotor_axis_alpha(alpha=axis_alpha,
pivot=array(self._param_pivot, float64), point=self.com)
# Pre-calculate the eigenframe rotation matrix.
two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen)
r22529 - in /branches/double_rotor: lib/frame_order/ specific_analyses/frame_order/ target_functions/