Author: bugman
Date: Thu Dec 8 13:39:49 2011
New Revision: 15059
URL: http://svn.gna.org/viewcvs/relax?rev=15059&view=rev
Log:
Added a script for testing the frame order analysis of the rotor model with
only PCSs.
Added:
branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order_pcs_only.py
- copied, changed from r15025,
branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order.py
Copied:
branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order_pcs_only.py
(from r15025,
branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order.py)
URL:
http://svn.gna.org/viewcvs/relax/branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order_pcs_only.py?p2=branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order_pcs_only.py&p1=branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order.py&r1=15025&r2=15059&rev=15059&view=diff
==============================================================================
---
branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order.py
(original)
+++
branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order_pcs_only.py
Thu Dec 8 13:39:49 2011
@@ -15,49 +15,73 @@
# Optimise.
self.optimisation()
- # The rotation matrix.
- R = zeros((3, 3), float64)
- euler_to_R_zyz(cdp.ave_pos_alpha, cdp.ave_pos_beta,
cdp.ave_pos_gamma, R)
- print("Rotation matrix:\n%s\n" % R)
- R = transpose(R)
- print("Inverted rotation:\n%s\n" % R)
-
- # The pivot point.
- pivot = array([ 37.254, 0.5, 16.7465])
-
# Load the original structure.
self.original_structure()
# Domain transformation.
- self.transform(R, pivot)
+ self.transform()
# Display in pymol.
- self.pymol_display(pivot)
+ self.pymol_display()
# Save the state.
- state.save('frame_order', force=True)
+ state.save('frame_order_pcs_only', force=True)
def optimisation(self):
"""Optimise the frame order model."""
- # The file paths.
- PATH_N_DOM = '..' + sep
- PATH_C_DOM = '.' + sep
-
# Create the data pipe.
pipe.create(pipe_name='frame order', pipe_type='frame order')
- # Load the tensors.
- script(PATH_N_DOM + 'tensors.py')
- script(PATH_C_DOM + 'tensors.py')
+ # Read the structures.
+ structure.read_pdb('1J7O_1st_NH.pdb', dir='..', set_mol_name='N-dom')
+ structure.read_pdb('1J7P_1st_NH_rot.pdb', dir='..',
set_mol_name='C-dom')
+
+ # Load the spins.
+ structure.load_spins('@N')
+ structure.load_spins('@H')
+
+ # Load the NH vectors.
+ structure.vectors(spin_id='@N', attached='H', ave=False)
+
+ # Set the values needed to calculate the dipolar constant.
+ value.set(1.041 * 1e-10, 'bond_length', spin_id="@N")
+ value.set('15N', 'heteronucleus', spin_id="@N")
+ value.set('1H', 'proton', spin_id="@N")
+
+ # Loop over the alignments.
+ ln = ['dy', 'tb', 'tm', 'er']
+ for i in range(len(ln)):
+ ## Load the RDCs.
+ #rdc.read(align_id=ln[i], file='rdc_%s.txt'%ln[i],
res_num_col=2, spin_name_col=5, data_col=6, error_col=7)
+
+ # The PCS.
+ pcs.read(align_id=ln[i], file='pcs_%s.txt'%ln[i], res_num_col=2,
spin_name_col=5, data_col=6, error_col=7)
+
+ # The temperature and field strength.
+ temperature(id=ln[i], temp=303)
+ frq.set(id=ln[i], frq=900e6)
+
+ # Load the N-domain tensors (the full tensors).
+ script('../tensors.py')
+
+ # Define the domains.
+ domain(id='N', spin_id=":1-78")
+ domain(id='C', spin_id=":80-144")
# The tensor domains and reductions.
full = ['Dy N-dom', 'Tb N-dom', 'Tm N-dom', 'Er N-dom']
red = ['Dy C-dom', 'Tb C-dom', 'Tm C-dom', 'Er C-dom']
for i in range(len(full)):
+ # Initalise the reduced tensor.
+ align_tensor.init(tensor=red[i], params=(0,0,0,0,0))
+
+ # Set the domain info.
align_tensor.set_domain(tensor=full[i], domain='N')
align_tensor.set_domain(tensor=red[i], domain='C')
+
+ # Specify which tensor is reduced.
align_tensor.reduction(full_tensor=full[i], red_tensor=red[i])
# Select the model.
@@ -66,17 +90,34 @@
# Set the reference domain.
frame_order.ref_domain('N')
+ # Set the initial pivot point.
+ pivot = array([ 37.254, 0.5, 16.7465])
+ frame_order.pivot(pivot)
+
+ # Set the paramagnetic centre.
+ paramag.centre(pos=[35.934, 12.194, -4.206])
+
+ # Check the minimum.
+ cdp.ave_pos_alpha = 1.2017352840543052
+ cdp.ave_pos_beta = 5.8477792871424867
+ cdp.ave_pos_gamma = 0.65969938507054027
+ cdp.axis_theta = 0.69832655838992175
+ cdp.axis_phi = 0.89069389677025046
+ cdp.cone_sigma_max = 0.52757207875029488
+ calc()
+ print cdp.chi2
+
# Optimise.
- grid_search(inc=3)
+ #grid_search(inc=5)
minimise('simplex', constraints=False)
- # Test Monte Carlo simulations.
- monte_carlo.setup(number=500)
- monte_carlo.create_data()
- monte_carlo.initial_values()
- minimise('simplex', constraints=False)
- eliminate()
- monte_carlo.error_analysis()
+ ## Test Monte Carlo simulations.
+ #monte_carlo.setup(number=500)
+ #monte_carlo.create_data()
+ #monte_carlo.initial_values()
+ #minimise('simplex', constraints=False)
+ #eliminate()
+ #monte_carlo.error_analysis()
def original_structure(self):
@@ -89,7 +130,7 @@
structure.read_pdb('1J7P_1st_NH.pdb', dir='..')
- def pymol_display(self, pivot):
+ def pymol_display(self):
"""Display the results in PyMOL."""
# Switch back to the main data pipe.
@@ -99,11 +140,8 @@
structure.read_pdb('1J7O_1st_NH.pdb', dir='..')
structure.read_pdb('1J7P_1st_NH_rot.pdb', dir='..')
- # Set the pivot point.
- frame_order.pivot(pivot)
-
# Create the cone PDB file.
- frame_order.cone_pdb(file='cone.pdb', force=True)
+ frame_order.cone_pdb(file='cone_pcs_only.pdb', force=True)
# Set the domains.
frame_order.domain_to_pdb(domain='N', pdb='1J7O_1st_NH.pdb')
@@ -112,11 +150,22 @@
# PyMOL.
pymol.view()
pymol.command('show spheres')
- pymol.cone_pdb('cone.pdb')
+ pymol.cone_pdb('cone_pcs_only.pdb')
- def transform(self, R, pivot):
+ def transform(self):
"""Transform the domain to the average position."""
+
+ # Switch back to the main data pipe.
+ pipe.switch('frame order')
+
+ # The rotation matrix.
+ R = zeros((3, 3), float64)
+ euler_to_R_zyz(cdp.ave_pos_alpha, cdp.ave_pos_beta,
cdp.ave_pos_gamma, R)
+ print("Rotation matrix:\n%s\n" % R)
+ R = transpose(R)
+ print("Inverted rotation:\n%s\n" % R)
+ pivot = cdp.pivot
# Create a special data pipe for the average rigid body position.
pipe.create(pipe_name='ave pos', pipe_type='frame order')
@@ -128,7 +177,7 @@
structure.rotate(R=R, origin=pivot)
# Write out the new PDB.
- structure.write_pdb('ave_pos', force=True)
+ structure.write_pdb('ave_pos_pcs_only', force=True)
# Execute the analysis.
r15059 - /branches/frame_order_testing/test_suite/shared_data/frame_order/rotor2/frame_order_pcs_only.py