Skip to content
Snippets Groups Projects
Commit 0fc54155 authored by Gonzalez, Miguel's avatar Gonzalez, Miguel Committed by Remi Perenon
Browse files

Debugging and modifying output variables

parent f496e99d
No related branches found
No related tags found
1 merge request!42Feature total structure factor
......@@ -101,26 +101,24 @@ class DynamicTotalStructureFactor(IJob):
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")
self._outputData.add("f(q,t)_inc_%s" % element,"surface", (nQShells,self._nFrames), axis="q|time", units="au")
self._outputData.add("s(q,f)_inc_%s" % element,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
self._outputData.add("f(q,t)_inc_weighted_%s" % element,"surface", (nQShells,self._nFrames), axis="q|time", units="au")
self._outputData.add("s(q,f)_inc_weighted_%s" % element,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
for pair in self._elementsPairs:
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_coh_weighted(q,t)_total","surface", (nQShells,self._nFrames), axis="q|time", units="au")
self._outputData.add("f_inc_weighted(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_coh_weighted(q,f)_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
self._outputData.add("s_inc_weighted(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")
self._outputData.add("f(q,t)_coh_%s%s" % pair,"surface", (nQShells,self._nFrames), axis="q|time" , units="au")
self._outputData.add("s(q,f)_coh_%s%s" % pair,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
self._outputData.add("f(q,t)_coh_weighted_%s%s" % pair,"surface", (nQShells,self._nFrames), axis="q|time" , units="au")
self._outputData.add("s(q,f)_coh_weighted_%s%s" % pair,"surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
self._outputData.add("f(q,t)_coh_total","surface", (nQShells,self._nFrames), axis="q|time", units="au")
self._outputData.add("f(q,t)_inc_total","surface", (nQShells,self._nFrames), axis="q|time", units="au")
self._outputData.add("f(q,t)_total","surface", (nQShells,self._nFrames), axis="q|time", units="au")
self._outputData.add("s(q,f)_coh_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
self._outputData.add("s(q,f)_inc_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
self._outputData.add("s(q,f)_total","surface", (nQShells,self._nOmegas), axis="q|omega", units="nm2/ps")
def run_step(self, index):
"""
......@@ -199,10 +197,10 @@ class DynamicTotalStructureFactor(IJob):
for pair in self._elementsPairs:
corr = correlation(rho[pair[0]],rho[pair[1]], average=1)
self._outputData["f_coh(q,t)_%s%s" % pair][index,:] += corr
self._outputData["f(q,t)_coh_%s%s" % pair][index,:] += corr
for k,v in disf.items():
self._outputData["f_inc(q,t)_%s" % k][index,:] += v
self._outputData["f(q,t)_inc_%s" % k][index,:] += v
def finalize(self):
"""
......@@ -214,44 +212,52 @@ class DynamicTotalStructureFactor(IJob):
nTotalAtoms = 0
for val in nAtomsPerElement.values():
nTotalAtoms += val
# Compute coherent functions and structure factor
for pair in self._elementsPairs:
b_coh1 = ELEMENTS[pair[0],"b_coherent"]
b_coh2 = ELEMENTS[pair[1],"b_coherent"]
bi = ELEMENTS[pair[0],"b_coherent"]
bj = ELEMENTS[pair[1],"b_coherent"]
ni = nAtomsPerElement[pair[0]]
nj = nAtomsPerElement[pair[1]]
ci = float(ni)/nTotalAtoms
cj = float(ni)/nTotalAtoms
cj = float(nj)/nTotalAtoms
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._outputData["f(q,t)_coh_%s%s" % pair][:] /= numpy.sqrt(ni*nj)
self._outputData["s(q,f)_coh_%s%s" % pair][:] = get_spectrum(self._outputData["f(q,t)_coh_%s%s" % pair],
self.configuration["instrument_resolution"]["time_window"],
self.configuration["instrument_resolution"]["time_step"],
axis=1)
self._outputData["f_coh_weighted(q,t)_total"][:] += self._outputData["f_coh(q,t)_%s%s" % pair][:] * numpy.sqrt(ci*cj) * b_coh1 * b_coh2
self._outputData["s_coh_weighted(q,f)_total"][:] += self._outputData["s_coh(q,f)_%s%s" % pair][:] * numpy.sqrt(ci*cj) * b_coh1 * b_coh2
self._outputData["f(q,t)_coh_weighted_%s%s" % pair][:] = self._outputData["f(q,t)_coh_%s%s" % pair][:] * numpy.sqrt(ci*cj) * bi * bj
self._outputData["s(q,f)_coh_weighted_%s%s" % pair][:] = self._outputData["s(q,f)_coh_%s%s" % pair][:] * numpy.sqrt(ci*cj) * bi * bj
self._outputData["f(q,t)_coh_total"][:] += self._outputData["f(q,t)_coh_weighted_%s%s" % pair][:]
self._outputData["s(q,f)_coh_total"][:] += self._outputData["s(q,f)_coh_weighted_%s%s" % pair][:]
# Compute incoherent functions and structure factor
for element, ni in nAtomsPerElement.items():
b_inc = ELEMENTS[element,"b_incoherent"]
bi = ELEMENTS[element,"b_incoherent2"]
ni = nAtomsPerElement[element]
ci = float(ni)/nTotalAtoms
self._outputData["f_inc(q,t)_%s" % element][:] /= ni
self._outputData["s_inc(q,f)_%s" % element][:] = get_spectrum(self._outputData["f_inc(q,t)_%s" % element],
self._outputData["f(q,t)_inc_%s" % element][:] /= ni
self._outputData["s(q,f)_inc_%s" % element][:] = get_spectrum(self._outputData["f(q,t)_inc_%s" % element],
self.configuration["instrument_resolution"]["time_window"],
self.configuration["instrument_resolution"]["time_step"],
axis=1)
self._outputData["f_inc_weighted(q,t)_total"][:] += self._outputData["f_inc(q,t)_%s" % element][:] * ci * b_inc**2
self._outputData["s_inc_weighted(q,f)_total"][:] += self._outputData["s_inc(q,f)_%s" % element][:] * ci * b_inc**2
self._outputData["f(q,t)_inc_weighted_%s" % element][:] = self._outputData["f(q,t)_inc_%s" % element][:] * ci * bi
self._outputData["s(q,f)_inc_weighted_%s" % element][:] = self._outputData["s(q,f)_inc_%s" % element][:] * ci * bi
self._outputData["f(q,t)_inc_total"][:] += self._outputData["f(q,t)_inc_weighted_%s" % element][:]
self._outputData["s(q,f)_inc_total"][:] += self._outputData["s(q,f)_inc_weighted_%s" % element][:]
# Compute total sum
self._outputData["f_tot(q,t)_total"][:] = self._outputData["f_coh_weighted(q,t)_total"][:] + self._outputData["f_inc_weighted(q,t)_total"][:]
self._outputData["s_tot(q,f)_total"][:] = self._outputData["s_coh_weighted(q,f)_total"][:] + self._outputData["s_inc_weighted(q,f)_total"][:]
# Compute total F(Q,t) = inc + coh
self._outputData["f(q,t)_total"][:] = self._outputData["f(q,t)_coh_total"][:] + self._outputData["f(q,t)_inc_total"][:]
self._outputData["s(q,f)_total"][:] = self._outputData["s(q,f)_coh_total"][:] + self._outputData["s(q,f)_inc_total"][:]
#
self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self._info)
self.configuration['trajectory']['instance'].close()
REGISTRY['ndtsf'] = NeutronDynamicTotalStructureFactor
\ No newline at end of file
REGISTRY['ndtsf'] = NeutronDynamicTotalStructureFactor
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment