specific_analyses.n_state

1 ############################################################################### 2 # # 3 # Copyright (C) 2007-2014 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """The N-state model or structural ensemble analysis user functions.""" 24 25 # Python module imports. 26 from math import acos, cos, pi 27 from numpy import array, dot, float64, zeros 28 from numpy.linalg import norm 29 30 # relax module imports. 31 from lib.errors import RelaxError 32 from lib.geometry.rotations import euler_to_R_zyz, two_vect_to_R 33 from lib.io import open_write_file 34 from lib.structure.cones import Iso_cone 35 from lib.structure.internal.object import Internal 36 from lib.structure.geometric import generate_vector_dist, generate_vector_residues 37 from lib.structure.represent.cone import cone_edge 38 from pipe_control.pipes import check_pipe 39 from pipe_control.structure.mass import centre_of_mass 40 41

42 -def CoM(pivot_point=None, centre=None):

43 """Centre of mass analysis. 44 45 This function does an analysis of the centre of mass (CoM) of the N states. This includes 46 calculating the order parameter associated with the pivot-CoM vector, and the associated 47 cone of motions. The pivot_point argument must be supplied. If centre is None, then the 48 CoM will be calculated from the selected parts of the loaded structure. Otherwise it will 49 be set to the centre arg. 50 51 @param pivot_point: The pivot point in the structural file(s). 52 @type pivot_point: list of float of length 3 53 @param centre: The optional centre of mass vector. 54 @type centre: list of float of length 3 55 """ 56 57 # Test if the current data pipe exists. 58 check_pipe() 59 60 # Set the pivot point. 61 cdp.pivot_point = pivot_point 62 63 # The centre has been supplied. 64 if centre: 65 cdp.CoM = centre 66 67 # Calculate from the structure file. 68 else: 69 cdp.CoM = centre_of_mass() 70 71 # Calculate the vector between the pivot and CoM points. 72 cdp.pivot_CoM = array(cdp.CoM, float64) - array(cdp.pivot_point, float64) 73 74 # Calculate the unit vector between the pivot and CoM points. 75 unit_vect = cdp.pivot_CoM / norm(cdp.pivot_CoM) 76 77 # Initilise some data structures. 78 R = zeros((3, 3), float64) 79 vectors = zeros((cdp.N, 3), float64) 80 81 # Loop over the N states. 82 for c in range(cdp.N): 83 # Generate the rotation matrix. 84 euler_to_R_zyz(cdp.alpha[c], cdp.beta[c], cdp.gamma[c], R) 85 86 # Rotate the unit vector. 87 vectors[c] = dot(R, unit_vect) 88 89 # Multiply by the probability. 90 vectors[c] = vectors[c] * cdp.probs[c] 91 92 # Average of the unit vectors. 93 cdp.ave_unit_pivot_CoM = sum(vectors) 94 95 # The length reduction. 96 cdp.ave_pivot_CoM_red = norm(cdp.ave_unit_pivot_CoM) 97 98 # The aveage pivot-CoM vector. 99 cdp.ave_pivot_CoM = norm(cdp.pivot_CoM) * cdp.ave_unit_pivot_CoM 100 101 # The full length rotated pivot-CoM vector. 102 cdp.full_ave_pivot_CoM = cdp.ave_pivot_CoM / cdp.ave_pivot_CoM_red 103 104 # The cone angle for diffusion on an axially symmetric cone. 105 cdp.theta_diff_on_cone = acos(cdp.ave_pivot_CoM_red) 106 cdp.S_diff_on_cone = (3.0*cos(cdp.theta_diff_on_cone)**2 - 1.0) / 2.0 107 108 # The cone angle and order parameter for diffusion in an axially symmetric cone. 109 cdp.theta_diff_in_cone = acos(2.*cdp.ave_pivot_CoM_red - 1.) 110 cdp.S_diff_in_cone = cos(cdp.theta_diff_in_cone) * (1 + cos(cdp.theta_diff_in_cone)) / 2.0 111 112 # Print out. 113 print("\n%-40s %-20s" % ("Pivot point:", repr(cdp.pivot_point))) 114 print("%-40s %-20s" % ("Moving domain CoM (prior to rotation):", repr(cdp.CoM))) 115 print("%-40s %-20s" % ("Pivot-CoM vector", repr(cdp.pivot_CoM))) 116 print("%-40s %-20s" % ("Pivot-CoM unit vector:", repr(unit_vect))) 117 print("%-40s %-20s" % ("Average of the unit pivot-CoM vectors:", repr(cdp.ave_unit_pivot_CoM))) 118 print("%-40s %-20s" % ("Average of the pivot-CoM vector:", repr(cdp.ave_pivot_CoM))) 119 print("%-40s %-20s" % ("Full length rotated pivot-CoM vector:", repr(cdp.full_ave_pivot_CoM))) 120 print("%-40s %-20s" % ("Length reduction from unity:", repr(cdp.ave_pivot_CoM_red))) 121 print("%-40s %.5f rad (%.5f deg)" % ("Cone angle (diffusion on a cone)", cdp.theta_diff_on_cone, cdp.theta_diff_on_cone / (2*pi) *360.)) 122 print("%-40s S_cone = %.5f (S^2 = %.5f)" % ("S_cone (diffusion on a cone)", cdp.S_diff_on_cone, cdp.S_diff_on_cone**2)) 123 print("%-40s %.5f rad (%.5f deg)" % ("Cone angle (diffusion in a cone)", cdp.theta_diff_in_cone, cdp.theta_diff_in_cone / (2*pi) *360.)) 124 print("%-40s S_cone = %.5f (S^2 = %.5f)" % ("S_cone (diffusion in a cone)", cdp.S_diff_in_cone, cdp.S_diff_in_cone**2)) 125 print("\n\n")

126 127

128 -def cone_pdb(cone_type=None, scale=1.0, file=None, dir=None, force=False):

129 """Create a PDB file containing a geometric object representing the various cone models. 130 131 Currently the only cone types supported are 'diff in cone' and 'diff on cone'. 132 133 134 @param cone_type: The type of cone model to represent. 135 @type cone_type: str 136 @param scale: The size of the geometric object is eqaul to the average pivot-CoM 137 vector length multiplied by this scaling factor. 138 @type scale: float 139 @param file: The name of the PDB file to create. 140 @type file: str 141 @param dir: The name of the directory to place the PDB file into. 142 @type dir: str 143 @param force: Flag which if set to True will cause any pre-existing file to be 144 overwritten. 145 @type force: int 146 """ 147 148 # Test if the cone models have been determined. 149 if cone_type == 'diff in cone': 150 if not hasattr(cdp, 'S_diff_in_cone'): 151 raise RelaxError("The diffusion in a cone model has not yet been determined.") 152 elif cone_type == 'diff on cone': 153 if not hasattr(cdp, 'S_diff_on_cone'): 154 raise RelaxError("The diffusion on a cone model has not yet been determined.") 155 else: 156 raise RelaxError("The cone type " + repr(cone_type) + " is unknown.") 157 158 # The number of increments for the filling of the cone objects. 159 inc = 20 160 161 # The rotation matrix. 162 R = zeros((3, 3), float64) 163 two_vect_to_R(array([0, 0, 1], float64), cdp.ave_pivot_CoM/norm(cdp.ave_pivot_CoM), R) 164 165 # The isotropic cone object. 166 if cone_type == 'diff in cone': 167 angle = cdp.theta_diff_in_cone 168 elif cone_type == 'diff on cone': 169 angle = cdp.theta_diff_on_cone 170 cone_obj = Iso_cone(angle) 171 172 # Create the structural object. 173 structure = Internal() 174 175 # Add a structure. 176 structure.add_molecule(name='cone') 177 178 # Alias the single molecule from the single model. 179 mol = structure.structural_data[0].mol[0] 180 181 # Add the pivot point. 182 mol.atom_add(pdb_record='HETATM', atom_num=1, atom_name='R', res_name='PIV', res_num=1, pos=cdp.pivot_point, element='C') 183 184 # Generate the average pivot-CoM vectors. 185 print("\nGenerating the average pivot-CoM vectors.") 186 sim_vectors = None 187 if hasattr(cdp, 'ave_pivot_CoM_sim'): 188 sim_vectors = cdp.ave_pivot_CoM_sim 189 res_num = generate_vector_residues(mol=mol, vector=cdp.ave_pivot_CoM, atom_name='Ave', res_name_vect='AVE', sim_vectors=sim_vectors, res_num=2, origin=cdp.pivot_point, scale=scale) 190 191 # Generate the cone outer edge. 192 print("\nGenerating the cone outer edge.") 193 cap_start_atom = mol.atom_num[-1]+1 194 cone_edge(mol=mol, cone_obj=cone_obj, res_name='CON', res_num=3, apex=cdp.pivot_point, R=R, scale=norm(cdp.pivot_CoM), inc=inc) 195 196 # Generate the cone cap, and stitch it to the cone edge. 197 if cone_type == 'diff in cone': 198 print("\nGenerating the cone cap.") 199 cone_start_atom = mol.atom_num[-1]+1 200 generate_vector_dist(mol=mol, res_name='CON', res_num=3, centre=cdp.pivot_point, R=R, phi_max_fn=cone_obj.phi_max, scale=norm(cdp.pivot_CoM), inc=inc) 201 202 # Create the PDB file. 203 print("\nGenerating the PDB file.") 204 pdb_file = open_write_file(file, dir, force=force) 205 structure.write_pdb(pdb_file) 206 pdb_file.close()

207

Source Code for Module specific_analyses.n_state_model.uf