DynamicCoherentStructureFactor.py 8.52 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 127 128 129
 
    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:
            traj = self.configuration['trajectory']['instance']
            
            qVectors = self.configuration["q_vectors"]["value"][shell]["q_vectors"]
130
                                                                        
131
            rho = {}
132
            for element in self.configuration['atom_selection']['unique_names']:
133 134 135 136 137 138 139
                rho[element] = numpy.zeros((self._nFrames, qVectors.shape[1]), dtype = numpy.complex64)

            # loop over the trajectory time steps
            for i, frame in enumerate(self.configuration['frames']['value']):
                
                conf = traj.configuration[frame]

140
                for element,idxs in self._indexesPerElement.items():
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155
                    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:
156
            corr = correlation(x[pair[0]],x[pair[1]], average=1)
157 158 159 160 161 162 163
            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 ...)
        """
        
164
        nAtomsPerElement = self.configuration['atom_selection'].get_natoms()
165
        for pair in self._elementsPairs:
166 167
            ni = nAtomsPerElement[pair[0]]
            nj = nAtomsPerElement[pair[1]]
168 169 170 171
            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"],
172
                                                                     axis=1)
173
        
174
        fqtTotal = weight(self.configuration["weights"].get_weights(),self._outputData,nAtomsPerElement,2,"f(q,t)_%s%s")
175 176 177

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