From 0bb0410eb4f84c3ea8974d8c2ac41edb805603c8 Mon Sep 17 00:00:00 2001 From: Eric Pellegrini Date: Wed, 26 Sep 2018 16:55:02 +0200 Subject: [PATCH] Implemented the dynamic total structure factor --- .../Jobs/DynamicCoherentStructureFactor.py | 7 +- .../Jobs/DynamicIncoherentStructureFactor.py | 7 +- .../Jobs/DynamicTotalStructureFactor.py | 249 ++++++++++++++++++ template_dtsf.py | 29 ++ 4 files changed, 286 insertions(+), 6 deletions(-) create mode 100644 Src/Framework/Jobs/DynamicTotalStructureFactor.py create mode 100755 template_dtsf.py diff --git a/Src/Framework/Jobs/DynamicCoherentStructureFactor.py b/Src/Framework/Jobs/DynamicCoherentStructureFactor.py index 7facccc..af2023f 100644 --- a/Src/Framework/Jobs/DynamicCoherentStructureFactor.py +++ b/Src/Framework/Jobs/DynamicCoherentStructureFactor.py @@ -138,9 +138,10 @@ class DynamicCoherentStructureFactor(IJob): #. x (any): The returned result(s) of run_step """ - for pair in self._elementsPairs: - corr = correlation(x[pair[0]],x[pair[1]], average=1) - self._outputData["f(q,t)_%s%s" % pair][index,:] += corr + if x is not None: + for pair in self._elementsPairs: + corr = correlation(x[pair[0]],x[pair[1]], average=1) + self._outputData["f(q,t)_%s%s" % pair][index,:] += corr def finalize(self): """ diff --git a/Src/Framework/Jobs/DynamicIncoherentStructureFactor.py b/Src/Framework/Jobs/DynamicIncoherentStructureFactor.py index 533130a..9ac7a78 100644 --- a/Src/Framework/Jobs/DynamicIncoherentStructureFactor.py +++ b/Src/Framework/Jobs/DynamicIncoherentStructureFactor.py @@ -132,9 +132,10 @@ class DynamicIncoherentStructureFactor(IJob): #. index (int): The index of the step.\n #. x (any): The returned result(s) of run_step """ - - for k,v in x.items(): - self._outputData["f(q,t)_%s" % k][index,:] += v + + if x is not None: + for k,v in x.items(): + self._outputData["f(q,t)_%s" % k][index,:] += v def finalize(self): """ diff --git a/Src/Framework/Jobs/DynamicTotalStructureFactor.py b/Src/Framework/Jobs/DynamicTotalStructureFactor.py new file mode 100644 index 0000000..0bb3f3e --- /dev/null +++ b/Src/Framework/Jobs/DynamicTotalStructureFactor.py @@ -0,0 +1,249 @@ +#MDANSE : Molecular Dynamics Analysis for Neutron Scattering Experiments +#------------------------------------------------------------------------------------------ +#Copyright (C) +#2015- Eric C. Pellegrini Institut Laue-Langevin +#BP 156 +#71 avenue des Martyrs +#38000 Grenoble Cedex 9 +#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 + +:author: Eric C. Pellegrini +''' + +import collections +import itertools + +import numpy + +from MDANSE import REGISTRY +from MDANSE import ELEMENTS +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 +from MDANSE.MolecularDynamics.Trajectory import read_atoms_trajectory + +class DynamicTotalStructureFactorError(Error): + pass + +class DynamicTotalStructureFactor(IJob): + """ + Computes the dynamic total structure factor for a set of atoms as the sum of the incoherent and coherent structure factors + """ + + label = "Dynamic Total Structure Factor" + + category = ('Analysis','Scattering',) + + ancestor = ["mmtk_trajectory"] + + 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'}}) + settings['atom_selection'] = ('atom_selection', {'dependencies':{'trajectory':'trajectory'}}) + 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'}}) + settings['output_files'] = ('output_files', {'formats':["netcdf","ascii"]}) + settings['running_mode'] = ('running_mode',{}) + + def initialize(self): + """ + Initialize the input parameters and analysis self variables + """ + + if not self.configuration['trajectory']['instance'].universe.is_periodic: + raise DynamicCoherentStructureFactorError('Cannot start %s analysis on non-periodic system' % self.label) + + if not self.configuration['q_vectors']['is_lattice']: + raise DynamicCoherentStructureFactorError('The Q vectors must be generated on a lattice to run %s analysis' % self.label) + + 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"] + + self._nOmegas = self._instrResolution['n_omegas'] + + self._outputData.add("q","line",self.configuration["q_vectors"]["shells"], units="inv_nm") + + self._outputData.add("time","line", self.configuration['frames']['duration'], units='ps') + self._outputData.add("time_window","line", self._instrResolution["time_window"], units="au") + + 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") + + self._elementsPairs = sorted(itertools.combinations_with_replacement(self.configuration['atom_selection']['unique_names'],2)) + self._indexesPerElement = self.configuration['atom_selection'].get_indexes() + + for element in self.configuration['atom_selection']['unique_names']: + self._outputData.add("f_inc(q,t)_%s" % element,"surface", (nQShells,self._nFrames), axis="q|time", units="au") + self._outputData.add("f_tot(q,t)_%s" % element,"surface", (nQShells,self._nFrames), axis="q|time", units="au") + self._outputData.add("s_inc(q,f)_%s" % element,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") + self._outputData.add("s_tot(q,f)_%s" % element,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") + + for pair in self._elementsPairs: + if pair[0] != pair[1]: + continue + self._outputData.add("f_coh(q,t)_%s%s" % pair,"surface", (nQShells,self._nFrames), axis="q|time" , units="au") + self._outputData.add("s_coh(q,f)_%s%s" % pair,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") + + self._outputData.add("f_coh(q,t)_total","surface", (nQShells,self._nFrames), axis="q|time", units="au") + self._outputData.add("f_inc(q,t)_total","surface", (nQShells,self._nFrames), axis="q|time", units="au") + self._outputData.add("f_tot(q,t)_total","surface", (nQShells,self._nFrames), axis="q|time", units="au") + + self._outputData.add("s_coh(q,f)_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") + self._outputData.add("s_inc(q,f)_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") + self._outputData.add("s_tot(q,f)_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps") + + 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'] + + nQVectors = self.configuration["q_vectors"]["value"][shell]["q_vectors"].shape[1] + + rho = {} + for element in self.configuration['atom_selection']['unique_names']: + rho[element] = numpy.zeros((self._nFrames, nQVectors), dtype = numpy.complex64) + + # loop over the trajectory time steps + for i, frame in enumerate(self.configuration['frames']['value']): + + qVectors = self.configuration["q_vectors"]["value"][shell]["q_vectors"] + + conf = traj.configuration[frame] + + for element,idxs in self._indexesPerElement.items(): + + selectedCoordinates = numpy.take(conf.array, idxs, axis=0) + rho[element][i,:] = numpy.sum(numpy.exp(1j*numpy.dot(selectedCoordinates, qVectors)),axis=0) + + + 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] + + 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]) + + temp = numpy.exp(1j*numpy.dot(series, qVectors)) + res = correlation(temp, axis=0, average=1) + + disf_per_q_shell[element] += res + + return index, (rho, disf_per_q_shell) + + + 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 + """ + + if x is not None: + rho,disf = x + + for pair in self._elementsPairs: + if pair[0] != pair[1]: + continue + corr = correlation(rho[pair[0]],rho[pair[1]], average=1) + self._outputData["f_coh(q,t)_%s%s" % pair][index,:] += corr + + for k,v in disf.items(): + self._outputData["f_inc(q,t)_%s" % k][index,:] += v + + def finalize(self): + """ + Finalizes the calculations (e.g. averaging the total term, output files creations ...) + """ + + nAtomsPerElement = self.configuration['atom_selection'].get_natoms() + for pair in self._elementsPairs: + if pair[0] != pair[1]: + continue + ni = nAtomsPerElement[pair[0]] + nj = nAtomsPerElement[pair[1]] + self._outputData["f_coh(q,t)_%s%s" % pair][:] /= numpy.sqrt(ni*nj) + self._outputData["s_coh(q,f)_%s%s" % pair][:] = get_spectrum(self._outputData["f_coh(q,t)_%s%s" % pair], + self.configuration["instrument_resolution"]["time_window"], + self.configuration["instrument_resolution"]["time_step"], + axis=1) + + for element, number in nAtomsPerElement.items(): + self._outputData["f_inc(q,t)_%s" % element][:] /= number + self._outputData["s_inc(q,f)_%s" % element][:] = get_spectrum(self._outputData["f_inc(q,t)_%s" % element], + self.configuration["instrument_resolution"]["time_window"], + self.configuration["instrument_resolution"]["time_step"], + axis=1) + + + for element, number in nAtomsPerElement.items(): + + b_coh = ELEMENTS[element,"b_coherent"] + b_inc = ELEMENTS[element,"b_incoherent"] + + self._outputData["f_tot(q,t)_%s" % element][:] += b_coh*self._outputData["f_coh(q,t)_%s%s" % (element,element)][:] + self._outputData["f_tot(q,t)_%s" % element][:] += b_inc*self._outputData["f_inc(q,t)_%s" % element][:] + + self._outputData["s_tot(q,f)_%s" % element][:] += b_coh*self._outputData["s_coh(q,f)_%s%s" % (element,element)][:] + self._outputData["s_tot(q,f)_%s" % element][:] += b_inc*self._outputData["s_inc(q,f)_%s" % element][:] + + + self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self._info) + + self.configuration['trajectory']['instance'].close() + +REGISTRY['dtsf'] = DynamicTotalStructureFactor + diff --git a/template_dtsf.py b/template_dtsf.py new file mode 100755 index 0000000..1a818ac --- /dev/null +++ b/template_dtsf.py @@ -0,0 +1,29 @@ +#!/usr/bin/python + +######################################################## +# This is an automatically generated MDANSE run script # +######################################################## + +from MDANSE import REGISTRY + +################################################################ +# Job parameters # +################################################################ + +parameters = {} +parameters['atom_selection'] = 'all' +parameters['atom_transmutation'] = None +parameters['frames'] = (0, 10, 1) +parameters['instrument_resolution'] = ('gaussian', {'mu': 0.0, 'sigma': 10.0}) +parameters['output_files'] = ('output', ['netcdf']) +parameters['q_vectors'] = ('spherical_lattice', {'width': 0.1, 'n_vectors': 50, 'shells': (0.1, 5, 0.1)}) +parameters['running_mode'] = ('monoprocessor', 1) +parameters['trajectory'] = '/home/pellegrini/git/mdanse/Data/Trajectories/MMTK/waterbox_in_periodic_universe.nc' +parameters['weights'] = 'b_coherent' + +################################################################ +# Setup and run the analysis # +################################################################ + +dtsf = REGISTRY['job']['dtsf']() +dtsf.run(parameters,status=True) -- 2.21.0