Author: bugman
Date: Mon Sep 22 15:08:06 2014
New Revision: 25956
URL: http://svn.gna.org/viewcvs/relax?rev=25956&view=rev
Log:
Created the specific_analyses.frame_order.geometric.generate_axis_system()
function.
This is now used by most parts of the frame order analysis to generate the
full 3D eigenframe of the
motions. Previously this was implemented each time the frame or major axis
was required. This
replicated and highly inconsistent code has been eliminated.
Modified:
branches/frame_order_cleanup/specific_analyses/frame_order/geometric.py
branches/frame_order_cleanup/specific_analyses/frame_order/uf.py
Modified:
branches/frame_order_cleanup/specific_analyses/frame_order/geometric.py
URL:
http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/specific_analyses/frame_order/geometric.py?rev=25956&r1=25955&r2=25956&view=diff
==============================================================================
--- branches/frame_order_cleanup/specific_analyses/frame_order/geometric.py
(original)
+++ branches/frame_order_cleanup/specific_analyses/frame_order/geometric.py
Mon Sep 22 15:08:06 2014
@@ -29,6 +29,7 @@
import sys
# relax module imports.
+from lib.errors import RelaxFault
from lib.frame_order.conversions 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
@@ -41,7 +42,7 @@
from lib.text.sectioning import subsection, subsubsection
from pipe_control.structure.mass import pipe_centre_of_mass
from specific_analyses.frame_order.data import domain_moving, generate_pivot
-from specific_analyses.frame_order.variables import MODEL_DOUBLE_ROTOR,
MODEL_FREE_ROTOR, MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR,
MODEL_ISO_CONE_TORSIONLESS, MODEL_LIST_FREE_ROTORS, MODEL_LIST_ISO_CONE,
MODEL_LIST_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE,
MODEL_PSEUDO_ELLIPSE_FREE_ROTOR, MODEL_PSEUDO_ELLIPSE_TORSIONLESS, MODEL_ROTOR
+from specific_analyses.frame_order.variables import MODEL_DOUBLE_ROTOR,
MODEL_FREE_ROTOR, MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR,
MODEL_ISO_CONE_TORSIONLESS, MODEL_LIST_DOUBLE, MODEL_LIST_FREE_ROTORS,
MODEL_LIST_ISO_CONE, MODEL_LIST_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE,
MODEL_PSEUDO_ELLIPSE_FREE_ROTOR, MODEL_PSEUDO_ELLIPSE_TORSIONLESS, MODEL_ROTOR
def add_axes(structure=None, representation=None, size=None, sims=False):
@@ -124,11 +125,7 @@
mol.atom_add(atom_num=1, atom_name='R', res_name='AXE',
res_num=1, pos=pivot1, element='C', pdb_record='HETATM')
# The axis system.
- axes = zeros((3, 3), float64)
- if sims:
- euler_to_R_zyz(cdp.eigen_alpha_sim[sim_indices[i]],
cdp.eigen_beta_sim[sim_indices[i]], cdp.eigen_gamma_sim[sim_indices[i]], axes)
- else:
- euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta,
cdp.eigen_gamma, axes)
+ axes = generate_axis_system(sim_index=sim_indices[i])
# Rotations and inversions.
axes = dot(T, axes)
@@ -184,24 +181,7 @@
pivot = generate_pivot(order=1, sim_index=sim_indices[i],
pdb_limit=True)
# The rotation matrix (rotation from the z-axis to the cone axis).
- R = zeros((3, 3), float64)
- if cdp.model in MODEL_LIST_PSEUDO_ELLIPSE:
- if sims:
- euler_to_R_zyz(cdp.eigen_alpha_sim[sim_indices[i]],
cdp.eigen_beta_sim[sim_indices[i]], cdp.eigen_gamma_sim[sim_indices[i]], R)
- else:
- euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta,
cdp.eigen_gamma, R)
- else:
- if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]:
- if sims:
- axis =
create_rotor_axis_alpha(alpha=cdp.axis_alpha_sim[sim_indices[i]],
pivot=pivot, point=com)
- else:
- axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha,
pivot=pivot, point=com)
- elif cdp.model in MODEL_LIST_ISO_CONE:
- if sims:
- axis =
create_rotor_axis_spherical(theta=cdp.axis_theta_sim[sim_indices[i]],
phi=cdp.axis_phi_sim[sim_indices[i]])
- else:
- axis = create_rotor_axis_spherical(theta=cdp.axis_theta,
phi=cdp.axis_phi)
- two_vect_to_R(array([0, 0, 1], float64), axis, R)
+ R = generate_axis_system(sim_index=sim_indices[i])
print(("Rotation matrix:\n%s" % R))
# The transformation.
@@ -358,27 +338,11 @@
# Alias the molecule.
mol = structure.get_molecule(mol_name, model=model_nums[i])
- # The single rotor models.
- if cdp.model not in [MODEL_DOUBLE_ROTOR]:
- # Generate the rotor axis.
- if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]:
- axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha,
pivot=pivot1, point=pipe_centre_of_mass(verbosity=0))
- elif cdp.model in MODEL_LIST_ISO_CONE:
- 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)
-
- # The double rotor models.
- else:
- # Generate the eigenframe of the motion.
- frame = zeros((3, 3), float64)
- euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma,
frame)
-
- # Add the z axis.
- axis = frame[:, 2]
+ # The axis system.
+ axes = generate_axis_system(sim_index=sim_indices[i])
# Transform the central axis.
- axis = dot(T, axis)
+ axis = dot(T, axes[:, 2])
# The label position.
pos = pivot1 + displacement*axis
@@ -450,21 +414,8 @@
com.append(pivot1)
# Generate the rotor axis.
- if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]:
- if sims:
-
axis.append(create_rotor_axis_alpha(alpha=cdp.axis_alpha_sim[sim_indices[i]],
pivot=pivot1, point=com[i]))
- else:
-
axis.append(create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot1,
point=com[i]))
- elif cdp.model in [MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR]:
- if sims:
-
axis.append(create_rotor_axis_spherical(theta=cdp.axis_theta_sim[sim_indices[i]],
phi=cdp.axis_phi_sim[sim_indices[i]]))
- else:
-
axis.append(create_rotor_axis_spherical(theta=cdp.axis_theta,
phi=cdp.axis_phi))
- else:
- if sims:
-
axis.append(create_rotor_axis_euler(alpha=cdp.eigen_alpha_sim[sim_indices[i]],
beta=cdp.eigen_beta_sim[sim_indices[i]],
gamma=cdp.eigen_gamma_sim[sim_indices[i]]))
- else:
-
axis.append(create_rotor_axis_euler(alpha=cdp.eigen_alpha,
beta=cdp.eigen_beta, gamma=cdp.eigen_gamma))
+ axes = generate_axis_system(sim_index=sim_indices[i])
+ axis.append(axes[:, 2])
# The size of the rotor, taking the 30 Angstrom cone
representation into account.
if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]:
@@ -505,11 +456,7 @@
com.append(pivot1)
# Generate the eigenframe of the motion.
- frame = zeros((3, 3), float64)
- if sims:
- euler_to_R_zyz(cdp.eigen_alpha_sim[sim_indices[i]],
cdp.eigen_beta_sim[sim_indices[i]], cdp.eigen_gamma_sim[sim_indices[i]],
frame)
- else:
- euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta,
cdp.eigen_gamma, frame)
+ frame = generate_axis_system(sim_index=sim_indices[i])
# Add the x and y axes.
axis.append(frame[:, 0])
@@ -765,3 +712,57 @@
pdb_file = open_write_file(file_root[i]+'.pdb', dir,
compress_type=compress_type, force=force)
structures[i].write_pdb(pdb_file)
pdb_file.close()
+
+
+def generate_axis_system(sim_index=None):
+ """Generate and return the full 3D axis system for the current frame
order model.
+
+ @keyword sim_index: The optional MC simulation index to obtain the frame
for a given simulation.
+ @type sim_index: None or int
+ @return: The full 3D axis system for the model.
+ @rtype: numpy rank-2, 3D float64 array
+ """
+
+ # Initialise.
+ axis = zeros(3, float64)
+ frame = zeros((3, 3), float64)
+
+ # The 1st pivot point.
+ pivot = generate_pivot(order=1, sim_index=sim_index, pdb_limit=True)
+
+ # The CoM of the system.
+ com = pipe_centre_of_mass(verbosity=0)
+
+ # The system for the rotor models.
+ if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]:
+ # Generate the axis.
+ if sim_index == None:
+ axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha,
pivot=pivot, point=com)
+ else:
+ axis =
create_rotor_axis_alpha(alpha=cdp.axis_alpha_sim[sim_index], pivot=pivot,
point=com)
+
+ # The system for the isotropic cones.
+ elif cdp.model in MODEL_LIST_ISO_CONE:
+ # Generate the axis.
+ if sim_index == None:
+ axis = create_rotor_axis_spherical(theta=cdp.axis_theta,
phi=cdp.axis_phi)
+ else:
+ axis =
create_rotor_axis_spherical(theta=cdp.axis_theta_sim[sim_index],
phi=cdp.axis_phi_sim[sim_index])
+
+ # Create a full normalised axis system.
+ two_vect_to_R(array([1, 0, 0], float64), axis, frame)
+
+ # The system for the pseudo-ellipses and double rotor.
+ elif cdp.model in MODEL_LIST_PSEUDO_ELLIPSE + MODEL_LIST_DOUBLE:
+ # Generate the eigenframe of the motion.
+ if sim_index == None:
+ euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma,
frame)
+ else:
+ euler_to_R_zyz(cdp.eigen_alpha_sim[sim_index],
cdp.eigen_beta_sim[sim_index], cdp.eigen_gamma_sim[sim_index], frame)
+
+ # Unknown model.
+ else:
+ raise RelaxFault
+
+ # Return the full eigenframe.
+ return frame
Modified: branches/frame_order_cleanup/specific_analyses/frame_order/uf.py
URL:
http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/specific_analyses/frame_order/uf.py?rev=25956&r1=25955&r2=25956&view=diff
==============================================================================
--- branches/frame_order_cleanup/specific_analyses/frame_order/uf.py
(original)
+++ branches/frame_order_cleanup/specific_analyses/frame_order/uf.py Mon
Sep 22 15:08:06 2014
@@ -40,7 +40,7 @@
from lib.warnings import RelaxWarning
from pipe_control import pipes
from specific_analyses.frame_order.checks import check_domain, check_model,
check_parameters, check_pivot
-from specific_analyses.frame_order.geometric import average_position,
create_ave_pos, create_geometric_rep
+from specific_analyses.frame_order.geometric import average_position,
create_ave_pos, create_geometric_rep, generate_axis_system
from specific_analyses.frame_order.optimisation import count_sobol_points
from specific_analyses.frame_order.parameters import assemble_param_vector,
update_model
from specific_analyses.frame_order.variables import MODEL_ISO_CONE,
MODEL_ISO_CONE_FREE_ROTOR, MODEL_ISO_CONE_TORSIONLESS, MODEL_LIST,
MODEL_LIST_FREE_ROTORS, MODEL_LIST_ISO_CONE, MODEL_LIST_PSEUDO_ELLIPSE,
MODEL_LIST_RESTRICTED_TORSION, MODEL_PSEUDO_ELLIPSE,
MODEL_PSEUDO_ELLIPSE_TORSIONLESS, MODEL_RIGID
@@ -125,37 +125,22 @@
angles = array([cdp.cone_theta_x, cdp.cone_theta_y, cone_sigma_max],
float64)
x, y, z = angles
- # The system for the isotropic cones.
- if cdp.model in MODEL_LIST_ISO_CONE:
- # Reconstruct the rotation axis.
- axis = zeros(3, float64)
- spherical_to_cartesian([1, cdp.axis_theta, cdp.axis_phi], axis)
-
- # Create a full normalised axis system.
- x_ax = array([1, 0, 0], float64)
- y_ax = cross(axis, x_ax)
- y_ax /= norm(y_ax)
- x_ax = cross(y_ax, axis)
- x_ax /= norm(x_ax)
- axes = transpose(array([x_ax, y_ax, axis], float64))
-
- # Start printout.
+ # The axis system.
+ axes = generate_axis_system()
+
+ # Start printout for the isotropic cones.
+ if cdp.model in MODEL_LIST_ISO_CONE:
print("\nOriginal parameters:")
print("%-20s %20.10f" % ("cone_theta", cdp.cone_theta))
print("%-20s %20.10f" % ("cone_sigma_max", cone_sigma_max))
print("%-20s %20.10f" % ("axis_theta", cdp.axis_theta))
print("%-20s %20.10f" % ("axis_phi", cdp.axis_phi))
- print("%-20s\n%s" % ("cone axis", axis))
+ print("%-20s\n%s" % ("cone axis", axes[:, 2]))
print("%-20s\n%s" % ("full axis system", axes))
print("\nPermutation '%s':" % permutation)
- # The system for the pseudo-ellipses.
- else:
- # Generate the eigenframe of the motion.
- frame = zeros((3, 3), float64)
- euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma,
frame)
-
- # Start printout.
+ # Start printout for the pseudo-ellipses.
+ else:
print("\nOriginal parameters:")
print("%-20s %20.10f" % ("cone_theta_x", cdp.cone_theta_x))
print("%-20s %20.10f" % ("cone_theta_y", cdp.cone_theta_y))
@@ -163,7 +148,7 @@
print("%-20s %20.10f" % ("eigen_alpha", cdp.eigen_alpha))
print("%-20s %20.10f" % ("eigen_beta", cdp.eigen_beta))
print("%-20s %20.10f" % ("eigen_gamma", cdp.eigen_gamma))
- print("%-20s\n%s" % ("eigenframe", frame))
+ print("%-20s\n%s" % ("eigenframe", axes))
print("\nPermutation '%s':" % permutation)
# The axis inversion structure.
@@ -239,10 +224,10 @@
# Permute the axes (pseudo-ellipses).
else:
- frame_new = transpose(array([inv[0]*frame[:, perm_axes[0]],
inv[1]*frame[:, perm_axes[1]], inv[2]*frame[:, perm_axes[2]]], float64))
+ axes_new = transpose(array([inv[0]*axes[:, perm_axes[0]],
inv[1]*axes[:, perm_axes[1]], inv[2]*axes[:, perm_axes[2]]], float64))
# Convert the permuted frame to Euler angles and store them.
- cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma =
R_to_euler_zyz(frame_new)
+ cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma =
R_to_euler_zyz(axes_new)
# End printout.
if cdp.model in MODEL_LIST_ISO_CONE:
@@ -262,7 +247,7 @@
print("%-20s %20.10f" % ("eigen_alpha", cdp.eigen_alpha))
print("%-20s %20.10f" % ("eigen_beta", cdp.eigen_beta))
print("%-20s %20.10f" % ("eigen_gamma", cdp.eigen_gamma))
- print("%-20s\n%s" % ("eigenframe", frame_new))
+ print("%-20s\n%s" % ("eigenframe", axes_new))
def pivot(pivot=None, order=1, fix=False):
r25956 - in /branches/frame_order_cleanup/specific_analyses/frame_order: geometric.py uf.py