Author: bugman Date: Wed Jun 4 09:04:34 2014 New Revision: 23627 URL: http://svn.gna.org/viewcvs/relax?rev=23627&view=rev Log: Updated the frame_order.py optimisation script for the small angle CaM rotor frame order test data. This now has the correct rotor torsion angle of 1 degree, and the spherical coordinates are now converted to the axis alpha parameter. Modified: branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor_small_angle/frame_order.py Modified: branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor_small_angle/frame_order.py URL: http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor_small_angle/frame_order.py?rev=23627&r1=23626&r2=23627&view=diff ============================================================================== --- branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor_small_angle/frame_order.py (original) +++ branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor_small_angle/frame_order.py Wed Jun 4 09:04:34 2014 @@ -1,7 +1,46 @@ -# Script for optimising the rotor frame order test model of CaM. +# Script for optimising the small angle (1 degree) rotor frame order test model of CaM. # Python module imports. -from numpy import array +from numpy import array, cross, float64, zeros +from numpy.linalg import norm + +# relax module imports. +from lib.frame_order.rotor_axis import create_rotor_axis_alpha +from lib.geometry.lines import closest_point_ax +from lib.geometry.coord_transform import spherical_to_cartesian +from lib.geometry.rotations import reverse_euler_zyz +from lib.geometry.vectors import vector_angle +from pipe_control.structure.mass import pipe_centre_of_mass + + +def alpha_angle(pivot=None, com=None, axis=None): + """Calculate and return the rotor alpha angle.""" + + # The CoM-pivot axis. + com_piv = com - pivot + com_piv /= norm(com_piv) + + # The mu_xy vector. + z_axis = array([0, 0, 1], float64) + mu_xy = cross(z_axis, com_piv) + mu_xy /= norm(mu_xy) + + # The alpha angle. + return vector_angle(mu_xy, axis, com_piv) + + +def shift_pivot(pivot_orig=None, com=None, axis=None): + """Shift the pivot to the closest point on the rotor axis to the CoM.)""" + + # The closest point. + pivot_new = closest_point_ax(line_pt=pivot_orig, axis=axis, point=com) + + # Printout. + print("\n%-20s%s" % ("Original pivot:", pivot_orig)) + print("%-20s%s" % ("New pivot:", pivot_new)) + + # Return the shifted pivot. + return pivot_new # The real parameter values. @@ -9,7 +48,12 @@ AVE_POS_ALPHA, AVE_POS_BETA, AVE_POS_GAMMA = [5.623468683852550, 0.435439748282942, 5.081265879629926] AXIS_THETA = 0.9600799785953431 AXIS_PHI = 4.0322755062196229 -CONE_SIGMA_MAX = 30.0 / 360.0 * 2.0 * pi +CONE_SIGMA_MAX = 1.0 / 360.0 * 2.0 * pi + +# Reconstruct the rotation axis. +AXIS = zeros(3, float64) +spherical_to_cartesian([1, AXIS_THETA, AXIS_PHI], AXIS) +print("Rotation axis: %s" % AXIS) # Create the data pipe. pipe.create(pipe_name='frame order', pipe_type='frame order') @@ -71,7 +115,7 @@ frame_order.ref_domain('N') # Set the initial pivot point. -pivot = array([ 37.254, 0.5, 16.7465]) +pivot = shift_pivot(pivot_orig=array([ 37.254, 0.5, 16.7465]), com=pipe_centre_of_mass(verbosity=0), axis=AXIS) frame_order.pivot(pivot, fix=True) # Set the paramagnetic centre. @@ -87,8 +131,7 @@ value.set(param='ave_pos_alpha', val=AVE_POS_ALPHA) value.set(param='ave_pos_beta', val=AVE_POS_BETA) value.set(param='ave_pos_gamma', val=AVE_POS_GAMMA) -value.set(param='axis_theta', val=AXIS_THETA) -value.set(param='axis_phi', val=AXIS_PHI) +value.set(param='axis_alpha', val=alpha_angle(pivot=pivot, com=pipe_centre_of_mass(verbosity=0), axis=AXIS)) value.set(param='cone_sigma_max', val=CONE_SIGMA_MAX) calc() @@ -105,6 +148,17 @@ # Optimise the pivot and model. frame_order.pivot(pivot, fix=False) minimise('simplex') + +# The distance from the optimised pivot and the rotation axis. +opt_piv = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z]) +print("\n\nOptimised pivot displacement: %s" % norm(pivot - opt_piv)) +pt = closest_point_ax(line_pt=pivot, axis=AXIS, point=opt_piv) +print("Distance from axis: %s\n" % norm(pt - opt_piv)) + +# Recreate the axis and compare to the original. +opt_axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=opt_piv, point=pipe_centre_of_mass(verbosity=0)) +print("Original axis: %s" % AXIS) +print("Optimised axis: %s" % opt_axis) # Test Monte Carlo simulations. monte_carlo.setup(number=5)