Author: bugman
Date: Wed Aug 4 11:26:58 2010
New Revision: 11401
URL: http://svn.gna.org/viewcvs/relax?rev=11401&view=rev
Log:
Added support for the isotropic cone model, that is with a torsion angle
restriction of sigma_max.
Modified:
1.3/maths_fns/frame_order.py
1.3/maths_fns/frame_order_matrix_ops.py
1.3/specific_fns/frame_order.py
Modified: 1.3/maths_fns/frame_order.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/maths_fns/frame_order.py?rev=11401&r1=11400&r2=11401&view=diff
==============================================================================
--- 1.3/maths_fns/frame_order.py (original)
+++ 1.3/maths_fns/frame_order.py Wed Aug 4 11:26:58 2010
@@ -31,7 +31,7 @@
from generic_fns.frame_order import print_frame_order_2nd_degree
from maths_fns.alignment_tensor import to_5D, to_tensor
from maths_fns.chi2 import chi2
-from maths_fns.frame_order_matrix_ops import
compile_2nd_matrix_iso_cone_free_rotor, compile_2nd_matrix_pseudo_ellipse,
compile_2nd_matrix_pseudo_ellipse_free_rotor,
compile_2nd_matrix_pseudo_ellipse_torsionless, reduce_alignment_tensor
+from maths_fns.frame_order_matrix_ops import compile_2nd_matrix_iso_cone,
compile_2nd_matrix_iso_cone_free_rotor, compile_2nd_matrix_pseudo_ellipse,
compile_2nd_matrix_pseudo_ellipse_free_rotor,
compile_2nd_matrix_pseudo_ellipse_torsionless, reduce_alignment_tensor
from maths_fns.rotation_matrix import euler_to_R_zyz as euler_to_R
from relax_errors import RelaxError
@@ -186,7 +186,7 @@
theta, phi, theta_cone = params
# Generate the 2nd degree Frame Order super matrix.
- self.frame_order_2nd =
compile_2nd_matrix_iso_cone(self.frame_order_2nd, self.rot, self.z_axis,
self.cone_axis, theta, phi, theta_cone)
+ self.frame_order_2nd =
compile_2nd_matrix_iso_cone_free_rotor(self.frame_order_2nd, self.rot,
self.z_axis, self.cone_axis, theta, phi, theta_cone)
# Make the Frame Order matrix contiguous.
self.frame_order_2nd = self.frame_order_2nd.copy()
@@ -208,8 +208,32 @@
return val
+ def func_iso_cone(self, params):
+ """Target function for isotropic cone model optimisation.
+
+ This function optimises against alignment tensors.
+
+ @param params: The vector of parameter values {beta, gamma, theta,
phi, s1} where the first 2 are the tensor rotation Euler angles, the next two
are the polar and azimuthal angles of the cone axis, and s1 is the isotropic
cone order parameter.
+ @type params: list of float
+ @return: The chi-squared or SSE value.
+ @rtype: float
+ """
+
+ # Unpack the parameters.
+ ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi,
cone_theta, sigma_max = params
+
+ # Generate the 2nd degree Frame Order super matrix.
+ frame_order_2nd = compile_2nd_matrix_iso_cone(self.frame_order_2nd,
self.rot, self.z_axis, self.cone_axis, axis_theta, axis_phi, cone_theta,
sigma_max)
+
+ # Reduce and rotate the tensors.
+ self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma,
frame_order_2nd)
+
+ # Return the chi-squared value.
+ return chi2(self.red_tensors, self.red_tensors_bc, self.red_errors)
+
+
def func_iso_cone_free_rotor(self, params):
- """Target function for isotropic cone model optimisation using the
alignment tensors.
+ """Target function for free rotor isotropic cone model optimisation.
This function optimises against alignment tensors.
@@ -233,7 +257,7 @@
def func_pseudo_ellipse(self, params):
- """Target function for pseudo-elliptic cone model optimisation using
alignment tensors.
+ """Target function for pseudo-elliptic cone model optimisation.
@param params: The vector of parameter values {alpha, beta, gamma,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y,
cone_sigma_max} where the first 3 are the average position rotation Euler
angles, the next 3 are the Euler angles defining the eigenframe, and the last
3 are the pseudo-elliptic cone geometric parameters.
@type params: list of float
@@ -255,7 +279,7 @@
def func_pseudo_ellipse_free_rotor(self, params):
- """Target function for free_rotor pseudo-elliptic cone model
optimisation using alignment tensors.
+ """Target function for free_rotor pseudo-elliptic cone model
optimisation.
@param params: The vector of parameter values {alpha, beta, gamma,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the
first 3 are the average position rotation Euler angles, the next 3 are the
Euler angles defining the eigenframe, and the last 2 are the free_rotor
pseudo-elliptic cone geometric parameters.
@type params: list of float
@@ -277,7 +301,7 @@
def func_pseudo_ellipse_torsionless(self, params):
- """Target function for torsionless pseudo-elliptic cone model
optimisation using alignment tensors.
+ """Target function for torsionless pseudo-elliptic cone model
optimisation.
@param params: The vector of parameter values {alpha, beta, gamma,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the
first 3 are the average position rotation Euler angles, the next 3 are the
Euler angles defining the eigenframe, and the last 2 are the torsionless
pseudo-elliptic cone geometric parameters.
@type params: list of float
@@ -299,7 +323,7 @@
def func_rigid(self, params):
- """Target function for rigid model optimisation using the alignment
tensors.
+ """Target function for rigid model optimisation.
This function optimises against alignment tensors. The Euler angles
for the tensor rotation are the 3 parameters optimised in this model.
Modified: 1.3/maths_fns/frame_order_matrix_ops.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/maths_fns/frame_order_matrix_ops.py?rev=11401&r1=11400&r2=11401&view=diff
==============================================================================
--- 1.3/maths_fns/frame_order_matrix_ops.py (original)
+++ 1.3/maths_fns/frame_order_matrix_ops.py Wed Aug 4 11:26:58 2010
@@ -59,6 +59,42 @@
matrix[2, 2] = fact * quad(part_int_daeg1_pseudo_ellipse_zz, -pi, pi,
args=(theta_x, theta_y, sigma_max), full_output=1)[0]
+def compile_2nd_matrix_iso_cone(matrix, R, z_axis, cone_axis, theta_axis,
phi_axis, cone_theta, sigma_max):
+ """Generate the rotated 2nd degree Frame Order matrix for the isotropic
cone.
+
+ The cone axis is assumed to be parallel to the z-axis in the eigenframe.
+
+ @param matrix: The Frame Order matrix, 2nd degree to be populated.
+ @type matrix: numpy 9D, rank-2 array
+ @param R: The rotation matrix to be populated.
+ @type R: numpy 3D, rank-2 array
+ @param z_axis: The molecular frame z-axis from which the cone axis
is rotated from.
+ @type z_axis: numpy 3D, rank-1 array
+ @param cone_axis: The storage structure for the cone axis.
+ @type cone_axis: numpy 3D, rank-1 array
+ @param theta_axis: The cone axis polar angle.
+ @type theta_axis: float
+ @param phi_axis: The cone axis azimuthal angle.
+ @type phi_axis: float
+ @param cone_theta: The cone opening angle.
+ @type cone_theta: float
+ @param sigma_max: The maximum torsion angle.
+ @type sigma_max: float
+ """
+
+ # Generate the cone axis from the spherical angles.
+ generate_vector(cone_axis, theta_axis, phi_axis)
+
+ # Populate the Frame Order matrix in the eigenframe.
+ populate_2nd_eigenframe_iso_cone(matrix, cone_theta, sigma_max)
+
+ # Average position rotation.
+ two_vect_to_R(z_axis, cone_axis, R)
+
+ # Rotate and return the frame order matrix.
+ return rotate_daeg(matrix, R)
+
+
def compile_2nd_matrix_iso_cone_free_rotor(matrix, R, z_axis, cone_axis,
theta_axis, phi_axis, s1):
"""Generate the rotated 2nd degree Frame Order matrix for the free rotor
isotropic cone.
@@ -1048,6 +1084,56 @@
matrix[2, 2] = (cos(angle) + 1.0) / 2.0
+def populate_2nd_eigenframe_iso_cone(matrix, tmax, smax):
+ """Populate the 2nd degree Frame Order matrix in the eigenframe for the
isotropic cone.
+
+ The cone axis is assumed to be parallel to the z-axis in the eigenframe.
+
+
+ @param matrix: The Frame Order matrix, 2nd degree.
+ @type matrix: numpy 9D, rank-2 array
+ @param tmax: The cone opening angle.
+ @type tmax: float
+ @param smax: The maximum torsion angle.
+ @type smax: float
+ """
+
+ # Zeros.
+ for i in range(9):
+ for j in range(9):
+ matrix[i, j] = 0.0
+
+ # The surface area normalisation factor.
+ fact = 1.0 / (smax * (1.0 - cos(theta_max)))
+
+ # Repetitive trig calculations.
+ sin_smax = sin(smax)
+ sin_2smax = sin(2.0*smax)
+ sin_tmax = sin(tmax)
+ sin_2tmax = sin(2.0*tmax)
+
+ # Diagonal.
+ matrix[0, 0] = fact * ((sin_2smax + 4.0*smax)*sin_2tmax +
8.0*sin_2smax*sin_tmax + (6.0*sin_2smax + 24.0*smax)*tmax)/64.0
+ matrix[1, 1] = fact * (sin_2smax*sin_2tmax + (8.0*sin_2smax +
32.0*smax)*sin_tmax + 6.0*sin_2smax*tmax)/64.0
+ matrix[2, 2] = fact * (sin_smax*sin_2tmax + 4.0*sin_smax*sin_tmax +
2.0*sin_smax*tmax)/8.0
+ matrix[3, 3] = matrix[1, 1]
+ matrix[4, 4] = matrix[0, 0]
+ matrix[5, 5] = matrix[2, 2]
+ matrix[6, 6] = matrix[2, 2]
+ matrix[7, 7] = matrix[2, 2]
+ matrix[8, 8] = fact * (smax*sin_2tmax + 2.0*smax*tmax)/4.0
+
+ # Off diagonal set 1.
+ matrix[0, 4] = matrix[4, 0] = fact * - ((sin_2smax -
4.0*smax)*sin_2tmax + 8.0*sin_2smax*sin_tmax + (6.0*sin_2smax -
24.0*smax)*tmax)/64.0
+ matrix[0, 8] = matrix[8, 0] = fact * - (smax*sin_2tmax -
2.0*smax*tmax)/8.0
+ matrix[4, 8] = matrix[8, 4] = matrix[0, 8]
+
+ # Off diagonal set 2.
+ matrix[1, 3] = matrix[3, 1] = fact * (sin_2smax*sin_2tmax +
(8.0*sin_2smax - 32.0*smax)*sin_tmax + 6.0*sin_2smax*tmax)/64.0
+ matrix[2, 6] = matrix[6, 2] = fact * (sin_smax*sin_2tmax -
2.0*sin_smax*tmax)/8.0
+ matrix[5, 7] = matrix[7, 5] = matrix[2, 6]
+
+
def populate_2nd_eigenframe_iso_cone_free_rotor(matrix, s1):
"""Populate the 2nd degree Frame Order matrix in the eigenframe for the
free rotor isotropic cone.
Modified: 1.3/specific_fns/frame_order.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/specific_fns/frame_order.py?rev=11401&r1=11400&r2=11401&view=diff
==============================================================================
--- 1.3/specific_fns/frame_order.py (original)
+++ 1.3/specific_fns/frame_order.py Wed Aug 4 11:26:58 2010
@@ -85,7 +85,7 @@
"""
# Initialise the parameter array using the tensor rotation Euler
angles (average domain position).
- if cdp.model in ['iso cone', 'iso cone, torsionless', 'iso cone,
free rotor']:
+ if cdp.model in ['iso cone, torsionless', 'iso cone, free rotor']:
param_vect = [cdp.ave_pos_beta, cdp.ave_pos_gamma]
else:
param_vect = [cdp.ave_pos_alpha, cdp.ave_pos_beta,
cdp.ave_pos_gamma]
@@ -106,8 +106,10 @@
param_vect.append(cdp.cone_theta_x)
param_vect.append(cdp.cone_theta_y)
- # Cone parameters - single isotropic order parameter.
- elif cdp.model in ['iso cone', 'iso cone, torsionless', 'iso cone,
free rotor']:
+ # Cone parameters - single isotropic angle or order parameter.
+ elif cdp.model in ['iso cone']:
+ param_vect.append(cdp.cone_theta)
+ elif cdp.model in ['iso cone, torsionless', 'iso cone, free rotor']:
param_vect.append(cdp.cone_s1)
# Cone parameters - torsion angle.
@@ -497,7 +499,7 @@
# Build the parameter name list.
if not len(cdp.params):
# The tensor rotation, or average domain position.
- if cdp.model not in ['iso cone', 'iso cone, torsionless', 'iso
cone, free rotor']:
+ if cdp.model not in ['iso cone, torsionless', 'iso cone, free
rotor']:
cdp.params.append('ave_pos_alpha')
cdp.params.append('ave_pos_beta')
cdp.params.append('ave_pos_gamma')
@@ -518,8 +520,10 @@
cdp.params.append('cone_theta_x')
cdp.params.append('cone_theta_y')
- # Cone parameters - single isotropic order parameter.
- elif cdp.model in ['iso cone', 'iso cone, torsionless', 'iso
cone, free rotor']:
+ # Cone parameters - single isotropic angle or order parameter.
+ elif cdp.model in ['iso cone']:
+ cdp.params.append('cone_angle')
+ elif cdp.model in ['iso cone, torsionless', 'iso cone, free
rotor']:
cdp.params.append('cone_s1')
# Cone parameters - torsion angle.
@@ -571,8 +575,7 @@
elif cdp.model in ['pseudo-ellipse, torsionless', 'pseudo-ellipse,
free rotor']:
ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha,
eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = param_vector
elif cdp.model == 'iso cone':
- ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1,
cone_sigma_max = param_vector
- ave_pos_alpha = 0.0
+ ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta,
axis_phi, cone_theta, cone_sigma_max = param_vector
elif cdp.model in ['iso cone, torsionless', 'iso cone, free rotor']:
ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 =
param_vector
ave_pos_alpha = 0.0
@@ -781,7 +784,7 @@
# Parameters.
if (set == 'all' or set == 'params') and hasattr(cdp, 'model'):
# Initialise the parameter array using the tensor rotation Euler
angles (average domain position).
- if cdp.model not in ['iso cone', 'iso cone, torsionless', 'iso
cone, free rotor']:
+ if cdp.model not in ['iso cone, torsionless', 'iso cone, free
rotor']:
names.append('ave_pos_alpha%s' % suffix)
names.append('ave_pos_beta%s' % suffix)
names.append('ave_pos_gamma%s' % suffix)
@@ -802,8 +805,10 @@
names.append('cone_theta_x%s' % suffix)
names.append('cone_theta_y%s' % suffix)
- # Cone parameters - single isotropic order parameter.
- elif cdp.model in ['iso cone', 'iso cone, torsionless', 'iso
cone, free rotor']:
+ # Cone parameters - single isotropic angle or order parameter.
+ elif cdp.model in ['iso cone']:
+ names.append('cone_theta%s' % suffix)
+ elif cdp.model in ['iso cone, torsionless', 'iso cone, free
rotor']:
names.append('cone_s1%s' % suffix)
# Cone parameters - torsion angle.
r11401 - in /1.3: maths_fns/frame_order.py maths_fns/frame_order_matrix_ops.py specific_fns/frame_order.py