This example demonatrates a routine procedure of calculating phonon density of states from an experimental NeXus data file for a powder vanadium sample measured at ARCS, a direct-geometry neutron chopper spectrometer at the Spallation Neutron Source (SNS), Oak Ridge National Lab.
Get python tools ready, this may take a while
In [ ]:
import os, numpy as np
import histogram.hdf as hh, histogram as H
from matplotlib import pyplot as plt
%matplotlib notebook
# %matplotlib inline
import mantid
from multiphonon.sqe import plot as plot_sqe
from multiphonon.ui.getdos import Context, NxsWizardStart, QEGridWizardStart, GetDOSWizStart
Create a context for getdos. It stores the processing parameters.
In [ ]:
context=Context()
Create a new working directory and change into it. All intermediate results and final outputs will be in this new directory.
Please modify the following path to suit your need!
In [ ]:
workdir = os.path.abspath('./V_Ei120meV')
!mkdir -p {workdir}
%cd {workdir}
If you want to reuse a previously-saved context, please uncomment the following cell and execute
In [ ]:
# context.from_yaml('getdos2-V_Ei120meV-context.yaml')
Build download command
In [ ]:
dest = 'ARCS_V_annulus.nxs'
url = "https://mcvine.ornl.gov/multiphonon/ARCS_V_annulus.nxs"
cmd = 'wget %r -O %r' % (url, dest)
print cmd
Download: this will take a while (can be a few minutes to an hour, depending on internet speed)
In [ ]:
%%time
!{cmd} >log.download 2>err.download
The following command should show the downloaded file "ARCS_V_annulus.nxs"
In [ ]:
ls
Set sample_nxs
In [ ]:
context.sample_nxs = './ARCS_V_annulus.nxs'
To start, we need to set the locations of the data files measured for the sample and empty can (for background correction), as well as the experimental conditions such as incident neutron energy (Ei, in meV) and sample temperature (T, in Kelvin). The following wizard help you go through these steps.
Please use these example inputs in the forms of the following wizard:
In [ ]:
NxsWizardStart(context).show()
Save configuration so you can reuse it
In [ ]:
context.to_yaml('./getdos2-V_Ei120meV-context.yaml')
Now we are ready to reduce the experimental data to obtain the dynamical structure factor, S(Q,E).
S(Q,E) spectra for both the sample and the empty can is the starting point for getdos processing.
Run the following wizard to define the E and Q axes so that S(Q,E) spectra can be obtained the INS raw data.
After you input Q-axis and E-axis parameters and continue, the data reduction will start to obtain I(Q,E) spectrum from the experimental data. This step takes a little while and the circle at the top-right corner of this notebook interface below the "logout" button will be filled, meaning the notebook is busy computing. Please be patient. Wait until it turns open to continue to the next step.
In [ ]:
QEGridWizardStart(context).show()
Parameters are saved in the work dir. Uncomment the script below to see.
In [ ]:
%%script bash
cat work/raw2iqe-sample.params
Plot sample IQE
In [ ]:
iqe = hh.load('work/iqe.h5')
plt.figure(figsize=(6,4))
plot_sqe(iqe)
# plt.xlim(0, 11)
plt.clim(0, 3e-3)
This is a plot of vanadium S(Q, E) histogram.
QEGridWizardStart(context).show()Now integreate over the Q (momentum transfer) axis to obtain energy spectrum I(E)
In [ ]:
iqe2 = iqe.copy()
I = iqe2.I; I[I!=I] = 0 # remove NaNs
IE = iqe2.sum('Q') # sum over Q
plt.figure(figsize=(6,4))
plt.plot(IE.energy, IE.I)
Save configuration so you can reuse it
In [ ]:
context.to_yaml('./getdos2-V_Ei120meV-context.yaml')
Phonon DOS will be obtained from SQE histogram by an iterative procedure where multiphonon and multiple scattering corrections are applied to the measured SQE spectrum, assuming incoherent approximation, and the corrected spectrum is then converted to DOS.
Please use the following inputs:
Click "Run" button.
In [ ]:
context.initdos = ''
GetDOSWizStart(context).show()
Save context
In [ ]:
context.to_yaml('./getdos2-V_Ei120meV-context.yaml')
Print context
In [ ]:
# print context
Results are saved in "work" directory
In [ ]:
ls work/
Plot the final result for DOS
In [ ]:
dos = hh.load('work/final-dos.h5')
plt.figure(figsize=(5,3))
plt.plot(dos.E, dos.I)
plt.xlabel('Energy (meV)')
# plt.xlim(0, 30)
plt.tight_layout()
More plotting utils are available
In [ ]:
from multiphonon.backward import plotutils as pu
In [ ]:
plt.figure(figsize=(5,3))
pu.plot_dos_iteration('work/')
In [ ]:
plt.figure(figsize=(6,4))
pu.plot_residual('work/')
In [ ]:
plt.figure(figsize=(8, 4))
pu.plot_intermediate_result_se('work/round-4')
In [ ]: