OrderParameter.py 9.59 KB
Newer Older
1
2
3
4
5
#MDANSE : Molecular Dynamics Analysis for Neutron Scattering Experiments
#------------------------------------------------------------------------------------------
#Copyright (C)
#2015- Eric C. Pellegrini Institut Laue-Langevin
#BP 156
6
7
#71 avenue des Martyrs
#38000 Grenoble Cedex 9
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#France
#pellegrini[at]ill.fr
#goret[at]ill.fr
#aoun[at]ill.fr
#
#This library is free software; you can redistribute it and/or
#modify it under the terms of the GNU Lesser General Public
#License as published by the Free Software Foundation; either
#version 2.1 of the License, or (at your option) any later version.
#
#This library is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#Lesser General Public License for more details.
#
#You should have received a copy of the GNU Lesser General Public
#License along with this library; if not, write to the Free Software
#Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

''' 
Created on Apr 10, 2015

30
:author: Bachir Aoun and Eric C. Pellegrini
31
32
33
'''

import collections
34
import os
35
36
37
38
39
40
41
42
43
44
45
46

import numpy

from Scientific.Geometry import Vector
from Scientific.Geometry.Transformation import angleFromSineAndCosine, Rotation 

from MDANSE.Mathematics.Signal import correlation
from MDANSE.Framework.Jobs.IJob import IJob
from MDANSE.MolecularDynamics.Trajectory import read_atoms_trajectory

class OrderParameter(IJob):
    """
47
48
49
50
    The combination of NMR and MD simulation data is very powerful in the study of conformational dynamics of proteins. NMR relaxation 
    spectroscopy is a unique approach for a site-specific investigation of both global tumbling and internal motions of proteins. 
    The molecular motions modulate the magnetic interactions  between the nuclear spins and lead for each nuclear spin to a relaxation 
    behaviour which reflects its environment. 
51
	
52
53
54
55
	The relationship between microscopic motions and measured spin relaxation rates is given by Redfield's theory. The relaxation 
	measurements probe the relaxation dynamics of a selected nuclear spin at only a few frequencies. Moreover, only a limited number 
	of independent observables are accessible. Hence, to relate relaxation data to protein dynamics,  a dynamical model for molecular 
	motions or a functional form, depending on a limited number of adjustable parameters, are required. 
56
    
57
58
59
    The generalized order parameter, indicates the degree of spatial restriction of the internal motions of a bond vector, while the 
    characteristic time is an effective correlation time, setting the time scale of the internal relaxation processes. The resulting 
    values range from 0 (completely disordered) to 1 (fully ordered). 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
    
    **Calculation:** \n
    angle at time T is calculated as the following: \n   
    .. math:: \\overrightarrow{vector} =  \\overrightarrow{direction} - \\overrightarrow{origin}
    .. math:: \phi(T = T_{1}-T_{0}) = arcos(  \\overrightarrow{vector(T_{1})} . \\overrightarrow{vector(T_{0})} )
    
    **Output:** \n      
    #. angular_correlation_legendre_1st: :math:`<cos(\phi(T))>`
    #. angular_correlation_legendre_2nd: :math:`<\\frac{1}{2}(3cos(\phi(T))^{2}-1)>`
    
    **Acknowledgement**\n
    AOUN Bachir, PELLEGRINI Eric
    """
    
    type = 'op'
    
    label = "Order parameter"

    category = ('Dynamics',)
    
    ancestor = "mmtk_trajectory"

82
    settings = collections.OrderedDict()
83
    settings['trajectory'] = ('mmtk_trajectory', {'default':os.path.join('..','..','..','Data','Trajectories', 'MMTK', 'protein_in_periodic_universe.nc')})
84
    settings['frames'] = ('frames', {'dependencies':{'trajectory':'trajectory'}})
85
86
87
    settings['axis_selection'] = ('atoms_list', {'dependencies':{'trajectory':'trajectory'},
                                                 'nAtoms':2,
                                                 'default':('C284H438N84O79S7',('C','C_beta'))})
88
89
90
91
    settings['reference_direction'] = ('vector', {'default':[0,0,1], 'notNull':True, 'normalize':True})
    settings['per_axis'] = ('boolean', {'label':"output contribution per axis", 'default':False})
    settings['output_files'] = ('output_files', {'formats':["netcdf","ascii"]})
    settings['running_mode'] = ('running_mode',{})
92
93
94
95
96
97
98
                
    def initialize(self):
        """
        Initialize the input parameters and analysis self variables
        """

        self._nFrames = self.configuration['frames']['number']
99
        self._nAxis = self.configuration['axis_selection']['n_values']
100
101
102
103
104
        
        self.numberOfSteps = self._nAxis

        self._outputData.add("times","line", self.configuration['frames']['time'], units='ps')

105
        self._outputData.add("axis_index","line", numpy.arange(self.configuration['axis_selection']['n_values']), units='au')
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
                        
        self._zAxis = Vector([0,0,1])
        refAxis = self.configuration['reference_direction']['value']
        axis = self._zAxis.cross(refAxis)

        theta = angleFromSineAndCosine(axis.length(),self._zAxis*refAxis)

        try:
            self._rotation = Rotation(axis,theta).tensor.array
        except ZeroDivisionError:
            self._doRotation = False
        else:
            self._doRotation = True

        self._outputData.add('p1',"line", (self._nFrames,), axis='times', units="au") 
        self._outputData.add('p2',"line", (self._nFrames,), axis='times', units="au") 
        self._outputData.add('s2',"line", (self._nAxis,), axis='times', units="au") 
            
        if self.configuration['per_axis']['value']:
            self._outputData.add('p1_per_axis',"surface", (self._nAxis,self._nFrames), axis='axis_index|times', units="au") 
            self._outputData.add('p2_per_axis',"surface", (self._nAxis,self._nFrames), axis='axis_index|times', units="au") 

    def run_step(self, index):
        """
130
        Runs a single step of the job.
131
 
132
133
134
135
        :param index: the index of the step.
        :type index: int
        :return: a 2-tuple whose 1st element is the index of the step and 2nd element resp. p1, p2, s2 vectors.
        :rtype: 2-tuple
136
137
        """

138
        e1, e2 = self.configuration['axis_selection']['atoms'][index]
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
        
        e1 = read_atoms_trajectory(self.configuration["trajectory"]["instance"],
                                   e1,
                                   first=self.configuration['frames']['first'],
                                   last=self.configuration['frames']['last']+1,
                                   step=self.configuration['frames']['step'])

        e2 = read_atoms_trajectory(self.configuration["trajectory"]["instance"],
                                   e2,
                                   first=self.configuration['frames']['first'],
                                   last=self.configuration['frames']['last']+1,
                                   step=self.configuration['frames']['step'])

        diff = e2 - e1

        modulus = numpy.sqrt(numpy.sum(diff**2,1))
        
        diff /= modulus[:,numpy.newaxis]
        
        # shape (3,n)
        tDiff = diff.T
        
        if self._doRotation:
            diff = numpy.dot(self._rotation,tDiff)
                
        costheta = numpy.dot(self._zAxis.array,tDiff)
        sintheta = numpy.sqrt(numpy.sum(numpy.cross(self._zAxis,tDiff,axisb=0)**2,1))
        
        cosphi = tDiff[0,:]/sintheta
        sinphi = tDiff[1,:]/sintheta

        tr2         = 3.0*costheta**2 - 1.
        cos2phi     = 2.0*cosphi**2-1.
        sin2phi     = 2.0*sinphi*cosphi
        cossintheta = costheta*sintheta
        sintheta_sq = sintheta**2
                
        # 1st order legendre polynomia
        p1 = correlation(costheta)

        # Formula for the 2nd legendre polynomia applied to spherical coordinates
        p2 = (0.25*correlation(tr2) + 
              3.00*correlation(cosphi*cossintheta) + 
              3.00*correlation(sinphi*cossintheta) + 
              0.75*correlation(cos2phi*sintheta_sq) + 
              0.75*correlation(sin2phi*sintheta_sq))

        # s2 calculation (s2 = lim (t->+inf) p2)
        s2 = (0.75 * (numpy.sum(cos2phi*sintheta_sq)**2 + numpy.sum(sin2phi*sintheta_sq)**2) + \
              3.00 * (numpy.sum(cosphi*cossintheta)**2 + numpy.sum(sinphi*cossintheta)**2) +
              0.25 *  numpy.sum(tr2)**2) / self._nFrames**2
        
        return index, (p1, p2, s2)

    def combine(self, index, x):
        """
195
196
197
198
199
200
201
        Combines/synchronizes the output of `run_step` method.
        
        :param index: the index of the step.
        :type index: int
        :param x: the output of `run_step` method.
        :type x: any python object
        """     
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
        
        p1, p2, s2 = x
        
        self._outputData['p1'] += p1
        self._outputData['p2'] += p2
        self._outputData['s2'][index] = s2
        
        if self.configuration['per_axis']['value']:
            self._outputData['p1_per_axis'][index,:] = p1
            self._outputData['p2_per_axis'][index,:] = p2
        
    def finalize(self):
        """
        Finalizes the calculations (e.g. averaging the total term, output files creations ...).
        """ 
             
        self._outputData['p1'] /= self._nAxis
        self._outputData['p2'] /= self._nAxis
                
221
        self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self._info)
222
223
        
        self.configuration['trajectory']['instance'].close()