DynamicIncoherentStructureFactor.py 7.91 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
# **************************************************************************
#
# MDANSE: Molecular Dynamics Analysis for Neutron Scattering Experiments
#
# @file      Src/Framework/Jobs/DynamicIncoherentStructureFactor.py
# @brief     Implements module/class/test DynamicIncoherentStructureFactor
#
# @homepage  https://mdanse.org
# @license   GNU General Public License v3 or higher (see LICENSE)
# @copyright Institut Laue Langevin 2013-now
# @authors   Scientific Computing Group at ILL (see AUTHORS)
#
# **************************************************************************

eric pellegrini's avatar
test  
eric pellegrini committed
15 16 17 18
import collections

import numpy

19
from MDANSE import REGISTRY
eric pellegrini's avatar
test  
eric pellegrini committed
20 21 22 23 24 25 26
from MDANSE.Framework.Jobs.IJob import IJob
from MDANSE.Mathematics.Arithmetic import weight
from MDANSE.Mathematics.Signal import correlation, get_spectrum
from MDANSE.MolecularDynamics.Trajectory import read_atoms_trajectory

class DynamicIncoherentStructureFactor(IJob):
    """
27 28
    Computes the dynamic incoherent structure factor S_inc(Q,w) for a set of atoms.
	It can be compared to experimental data e.g. the quasielastic scattering due to diffusion processes.
eric pellegrini's avatar
test  
eric pellegrini committed
29 30 31 32
    """
    
    label = "Dynamic Incoherent Structure Factor"

33
    category = ('Analysis','Scattering',)
eric pellegrini's avatar
test  
eric pellegrini committed
34
    
35
    ancestor = ["mmtk_trajectory","molecular_viewer"]
eric pellegrini's avatar
test  
eric pellegrini committed
36
    
37 38 39 40 41
    settings = collections.OrderedDict()
    settings['trajectory']=('mmtk_trajectory',{})
    settings['frames']=('frames', {"dependencies":{'trajectory':'trajectory'}})
    settings['instrument_resolution'] = ('instrument_resolution',{"dependencies":{'trajectory':'trajectory', 'frames' : 'frames'}})
    settings['q_vectors'] = ('q_vectors',{"dependencies":{'trajectory':'trajectory'}})
42 43 44
    settings['atom_selection']=('atom_selection',{"dependencies":{'trajectory':'trajectory'}})
    settings['grouping_level']=('grouping_level',{"dependencies":{'trajectory':'trajectory','atom_selection':'atom_selection', 'atom_transmutation':'atom_transmutation'}})
    settings['atom_transmutation']=('atom_transmutation',{"dependencies":{'trajectory':'trajectory', 'atom_selection':'atom_selection'}})
45
    settings['projection']=('projection', {"label":"project coordinates"})
46
    settings['weights'] = ('weights', {'default':'b_incoherent2',"dependencies":{'trajectory':'trajectory','atom_selection':'atom_selection', 'atom_transmutation':'atom_transmutation'}})
47 48
    settings['output_files']=('output_files', {"formats":["netcdf","ascii"]})
    settings['running_mode']=('running_mode',{})
49
                    
eric pellegrini's avatar
test  
eric pellegrini committed
50 51 52 53 54
    def initialize(self):
        """
        Initialize the input parameters and analysis self variables
        """

55
        self.numberOfSteps = self.configuration["q_vectors"]["n_shells"]
eric pellegrini's avatar
test  
eric pellegrini committed
56 57 58 59 60 61 62

        self._nQShells = self.configuration["q_vectors"]["n_shells"]

        self._nFrames = self.configuration['frames']['number']
        
        self._instrResolution = self.configuration["instrument_resolution"]
        
63
        self._nOmegas = self._instrResolution['n_omegas']
eric pellegrini's avatar
test  
eric pellegrini committed
64

65
        self._outputData.add("q","line", self.configuration["q_vectors"]["shells"], units="inv_nm") 
eric pellegrini's avatar
test  
eric pellegrini committed
66

67
        self._outputData.add("time","line", self.configuration['frames']['duration'], units='ps')
68
        self._outputData.add("time_window","line", self._instrResolution["time_window"], axis="time", units="au") 
eric pellegrini's avatar
test  
eric pellegrini committed
69

70 71
        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") 
eric pellegrini's avatar
test  
eric pellegrini committed
72
                        
73
        for element in self.configuration['atom_selection']['unique_names']:
74
            self._outputData.add("f(q,t)_%s" % element,"surface", (self._nQShells,self._nFrames)     , axis="q|time", units="au")                                                 
75
            self._outputData.add("s(q,f)_%s" % element,"surface", (self._nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") 
eric pellegrini's avatar
test  
eric pellegrini committed
76

77
        self._outputData.add("f(q,t)_total","surface", (self._nQShells,self._nFrames)     , axis="q|time", units="au")                                                 
78
        self._outputData.add("s(q,f)_total","surface", (self._nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") 
eric pellegrini's avatar
test  
eric pellegrini committed
79 80 81 82 83 84 85 86 87 88 89 90
    
    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. 
            #. atomicSF (numpy.array): The atomic structure factor
        """
        
91
        shell = self.configuration["q_vectors"]["shells"][index]
eric pellegrini's avatar
test  
eric pellegrini committed
92
        
93 94 95 96
        if not shell in self.configuration["q_vectors"]["value"]:
            return index, None
            
        else:
eric pellegrini's avatar
test  
eric pellegrini committed
97

98
            qVectors = self.configuration["q_vectors"]["value"][shell]["q_vectors"]
eric pellegrini's avatar
test  
eric pellegrini committed
99

100 101 102 103 104 105 106 107 108 109

            disf_per_q_shell = {}
            for element in self.configuration['atom_selection']['unique_names']:
                disf_per_q_shell[element] = numpy.zeros((self._nFrames,), dtype = numpy.float)

            for i,atom_indexes in enumerate(self.configuration['atom_selection']["indexes"]):

                masses = self.configuration['atom_selection']["masses"][i]

                element = self.configuration['atom_selection']["names"][i]
eric pellegrini's avatar
test  
eric pellegrini committed
110
                        
111 112 113 114 115 116 117 118 119 120 121 122 123
                series = read_atoms_trajectory(self.configuration["trajectory"]["instance"],
                                               atom_indexes,
                                               first=self.configuration['frames']['first'],
                                               last=self.configuration['frames']['last']+1,
                                               step=self.configuration['frames']['step'],
                                               weights=[masses])
                
                series = self.configuration['projection']["projector"](series)
                                                                                    
                rho = numpy.exp(1j*numpy.dot(series, qVectors))
                res = correlation(rho, axis=0, average=1)
                    
                disf_per_q_shell[element] += res
eric pellegrini's avatar
test  
eric pellegrini committed
124
        
125
        return index, disf_per_q_shell
eric pellegrini's avatar
test  
eric pellegrini committed
126 127 128 129 130 131 132 133 134 135
    
    
    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
        """
                
136 137 138
        for k,v in x.items():
            self._outputData["f(q,t)_%s" % k][index,:] += v
                        
eric pellegrini's avatar
test  
eric pellegrini committed
139 140 141 142 143
    def finalize(self):
        """
        Finalizes the calculations (e.g. averaging the total term, output files creations ...)
        """
        
144 145
        nAtomsPerElement = self.configuration['atom_selection'].get_natoms()
        for element, number in nAtomsPerElement.items():
eric pellegrini's avatar
test  
eric pellegrini committed
146 147 148 149 150 151
            self._outputData["f(q,t)_%s" % element][:] /= number
            self._outputData["s(q,f)_%s" % element][:] = get_spectrum(self._outputData["f(q,t)_%s" % element],
                                                                      self.configuration["instrument_resolution"]["time_window"],
                                                                      self.configuration["instrument_resolution"]["time_step"],
                                                                      axis=1)

152 153 154
        weights = self.configuration["weights"].get_weights()
                
        self._outputData["f(q,t)_total"][:] = weight(weights,self._outputData,nAtomsPerElement,1,"f(q,t)_%s")
eric pellegrini's avatar
test  
eric pellegrini committed
155
        
156
        self._outputData["s(q,f)_total"][:] = weight(weights,self._outputData,nAtomsPerElement,1,"s(q,f)_%s")
eric pellegrini's avatar
test  
eric pellegrini committed
157
    
158
        self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self._info)
eric pellegrini's avatar
test  
eric pellegrini committed
159
        
160 161 162
        self.configuration['trajectory']['instance'].close()
        
REGISTRY['disf'] = DynamicIncoherentStructureFactor