DynamicCoherentStructureFactor.py 8.65 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: Eric C. Pellegrini
31 32 33 34 35 36 37
'''

import collections
import itertools

import numpy

38
from MDANSE import REGISTRY
39 40 41 42 43 44 45 46 47 48
from MDANSE.Core.Error import Error
from MDANSE.Framework.Jobs.IJob import IJob
from MDANSE.Mathematics.Arithmetic import weight
from MDANSE.Mathematics.Signal import correlation, get_spectrum

class DynamicCoherentStructureFactorError(Error):
    pass

class DynamicCoherentStructureFactor(IJob):
    """
49 50 51
    Computes the dynamic coherent structure factor S_coh(Q,w) for a set of atoms.
	It can be compared to experimental data e.g. the energy-integrated, static structure factor S_coh(Q)
	or the dispersion and intensity of phonons.
52 53 54 55
    """
    
    label = "Dynamic Coherent Structure Factor"

56
    category = ('Analysis','Scattering',)
57
    
58
    ancestor = ["mmtk_trajectory","molecular_viewer"]
59

60 61 62
    settings = collections.OrderedDict()
    settings['trajectory'] = ('mmtk_trajectory',{})
    settings['frames'] = ('frames', {'dependencies':{'trajectory':'trajectory'}})
63
    settings['instrument_resolution'] = ('instrument_resolution',{'dependencies':{'trajectory':'trajectory','frames' : 'frames'}})
64 65
    settings['q_vectors'] = ('q_vectors',{'dependencies':{'trajectory':'trajectory'}})
    settings['atom_selection'] = ('atom_selection', {'dependencies':{'trajectory':'trajectory'}})
66 67
    settings['atom_transmutation'] = ('atom_transmutation', {'dependencies':{'trajectory':'trajectory','atom_selection':'atom_selection'}})
    settings['weights'] = ('weights', {'default':'b_coherent',"dependencies":{'trajectory':'trajectory','atom_selection':'atom_selection', 'atom_transmutation':'atom_transmutation'}})
68 69
    settings['output_files'] = ('output_files', {'formats':["netcdf","ascii"]})
    settings['running_mode'] = ('running_mode',{})
70 71 72 73 74 75 76
    
    def initialize(self):
        """
        Initialize the input parameters and analysis self variables
        """

        if not self.configuration['trajectory']['instance'].universe.is_periodic:
77
            raise DynamicCoherentStructureFactorError('Cannot start %s analysis on non-periodic system' % self.label)
78 79
        
        if not self.configuration['q_vectors']['is_lattice']:
80
            raise DynamicCoherentStructureFactorError('The Q vectors must be generated on a lattice to run %s analysis' % self.label)
81 82 83 84 85 86 87 88 89
        
        self.numberOfSteps = self.configuration['q_vectors']['n_shells']
        
        nQShells = self.configuration["q_vectors"]["n_shells"]

        self._nFrames = self.configuration['frames']['number']
        
        self._instrResolution = self.configuration["instrument_resolution"]
        
90
        self._nOmegas = self._instrResolution['n_omegas']
91 92 93

        self._outputData.add("q","line",self.configuration["q_vectors"]["shells"], units="inv_nm") 

94
        self._outputData.add("time","line", self.configuration['frames']['duration'], units='ps')
95
        self._outputData.add("time_window","line", self._instrResolution["time_window"], units="au") 
96

97 98
        self._outputData.add("omega","line", self._instrResolution["omega"], units='rad/ps')
        self._outputData.add("omega_window","line", self._instrResolution["omega_window"], axis="omega", units="au") 
99
                                
100 101
        self._elementsPairs = sorted(itertools.combinations_with_replacement(self.configuration['atom_selection']['unique_names'],2))
        self._indexesPerElement = self.configuration['atom_selection'].get_indexes()
102 103

        for pair in self._elementsPairs:
104 105
            self._outputData.add("f(q,t)_%s%s" % pair,"surface", (nQShells,self._nFrames),axis="q|time", units="au")                                                 
            self._outputData.add("s(q,f)_%s%s" % pair,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") 
106 107

        self._outputData.add("f(q,t)_total","surface", (nQShells,self._nFrames), axis="q|time", units="au")                                                 
108
        self._outputData.add("s(q,f)_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") 
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
 
    def run_step(self, index):
        """
        Runs a single step of the job.\n
 
        :Parameters:
            #. index (int): The index of the step.
        :Returns:
            #. index (int): The index of the step. 
            #. rho (numpy.array): The exponential part of I(k,t)
        """
        
        shell = self.configuration["q_vectors"]["shells"][index]
        
        if not shell in self.configuration["q_vectors"]["value"]:
            return index, None
            
        else:
127

128 129
            traj = self.configuration['trajectory']['instance']
            
130 131
            nQVectors = self.configuration["q_vectors"]["value"][shell]["q_vectors"].shape[1]
                                                                                    
132
            rho = {}
133
            for element in self.configuration['atom_selection']['unique_names']:
134
                rho[element] = numpy.zeros((self._nFrames, nQVectors), dtype = numpy.complex64)
135 136 137

            # loop over the trajectory time steps
            for i, frame in enumerate(self.configuration['frames']['value']):
138 139 140
            
                qVectors = self.configuration["q_vectors"]["value"][shell]["q_vectors"]
                                
141 142
                conf = traj.configuration[frame]

143
                for element,idxs in self._indexesPerElement.items():
144

145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
                    selectedCoordinates = numpy.take(conf.array, idxs, axis=0)
                    rho[element][i,:] = numpy.sum(numpy.exp(1j*numpy.dot(selectedCoordinates, qVectors)),axis=0)

            return index, rho
    
    
    def combine(self, index, x):
        """
        Combines returned results of run_step.\n
        :Parameters:
            #. index (int): The index of the step.\n
            #. x (any): The returned result(s) of run_step
        """
               
        for pair in self._elementsPairs:
160
            corr = correlation(x[pair[0]],x[pair[1]], average=1)
161 162 163 164 165 166 167
            self._outputData["f(q,t)_%s%s" % pair][index,:] += corr
            
    def finalize(self):
        """
        Finalizes the calculations (e.g. averaging the total term, output files creations ...)
        """
        
168
        nAtomsPerElement = self.configuration['atom_selection'].get_natoms()
169
        for pair in self._elementsPairs:
170 171
            ni = nAtomsPerElement[pair[0]]
            nj = nAtomsPerElement[pair[1]]
172 173 174 175
            self._outputData["f(q,t)_%s%s" % pair][:] /= numpy.sqrt(ni*nj)
            self._outputData["s(q,f)_%s%s" % pair][:] = get_spectrum(self._outputData["f(q,t)_%s%s" % pair],
                                                                     self.configuration["instrument_resolution"]["time_window"],
                                                                     self.configuration["instrument_resolution"]["time_step"],
176
                                                                     axis=1)
177
        
178
        fqtTotal = weight(self.configuration["weights"].get_weights(),self._outputData,nAtomsPerElement,2,"f(q,t)_%s%s")
179 180 181

        self._outputData["f(q,t)_total"][:] = fqtTotal
        
182
        sqfTotal = weight(self.configuration["weights"].get_weights(),self._outputData,nAtomsPerElement,2,"s(q,f)_%s%s")
183 184
        self._outputData["s(q,f)_total"][:] = sqfTotal
    
185
        self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self._info)
186 187 188
        
        self.configuration['trajectory']['instance'].close()     
  
189
REGISTRY['dcsf'] = DynamicCoherentStructureFactor
190