>> model % display model informationTo create a model, use the ifitmakefunc dialogue window (see below), or instantiate an iFunc object.
>> disp(model) % idem, extensive information
>> plot(model) % plot the model with its default settings
>> model(p, x, y ...) % evaluate the model with parameters p, and axes x,y,...
>> model([], x, y ...) % evaluate the model with axes x,y,... and automatic parameter guess
>> model('guess', x, y ...) % idem
>> a(model, p) % evaluate the model onto the iData object 'a' axes with parameters p
>> fits(a, model, p) % fit the model onto the iData object 'a'
>> fits(model, a, p) % same as above
>> save(model, 'filename', 'mat') % save the Model as a Matlab file. orher formats are possible (YAML, m, JSON, XML...)
>> iFunc('filename') % import a Model stored as m, YAML, MAT, JSON file format.
>> model = gauss;You can edit their code to see how to define new models:
>> model = gauss+lorz;
>> model = iFunc('p(1)*x+p(2)')
>> model = iFunc('a=p(1); b=p(2); signal=a*x+b')
>> edit gauss % edit the function definition (from a file)
>> edit voigt
>> edit(voig) % edit the object definition
Function 
Description 
Dimensionality 
Parameters 
allometric 
Allometric (power/asymptotic
law) 
1D 
Amplitude Offset Exponent
BackGround 
bigauss 
Asymmetric Gaussian 
1D  Amplitude Centre HalfWidth1 HalfWidth2 Background 
bilorz 
Asymmetric Lorentzian 
1D  Amplitude Centre HalfWidth1 HalfWidth2 Background 
bose  Bose factor 
1D  Tau [h/2pi/kT] in "x" units 
dho 
Damped harmonic oscillator 
1D  Amplitude Centre HalfWidth
Background Temperature (in
"x" unit) 
dirac 
Dirac peak 
1D 
Amplitude Centre 
doseresp 
Doseresponse curve with
variable Hill slope. This is a sigmoid or Sshaped. 
1D 
Amplitude Center Slope
BackGround 
expon 
Exponential
decay 
1D  Amplitude Tau Background 
expstretched 
Exponential
 Stretched 
1D  Amplitude Tau Exponent
Background 
gauss  Gaussian

1D  Amplitude Centre HalfWidth
Background 
green 
Green function 
1D 
Amplitude Centre HalfWidth Background 
heaviside 
Heaviside (gap) The GapSide indicates raising (+) or falling () gap. 
1D 
Amplitude Centre GapSide
Background 
langevin 
Langevin function for
magnetic polarization 
1D 
Amplitude Center Width
BackGround 
laplace 
Laplace
distribution function 
1D  Amplitude Center Width BackGround 
lognormal 
LogNormal
distribution 
1D  Amplitude Center Width
BackGround 
lorz 
Lorentzian
(aka Cauchy) 
1D  Amplitude Centre HalfWidth Background 
ngauss 
multiple Gaussian 
n*1D 

nlorz 
multiple Lorentzian 
n*1D 

pareto 
Pareto distribution function 
1D 
Amplitude Exponent Width
BackGround 
poisson 
Poisson distribution WARNING: The 'x' axis is assumed to be an integer array (counts) 
1D 
Amplitude Center BackGround 
pseudovoigt 
Pseudo Voigt 
1D 
Amplitude Center Width
BackGround LorentzianRatio 
quadline 
Quadratic
line (parabola) 
1D  Quadratic Linear Constant 
sigmoid 
Sigmoidal
Sshaped curve (similar to Dose Response) 
1D  Amplitude Center Width BackGround 
sine 
Sine
function 
1D  Amplitude Phase_Shift Period
BackGround 
sinedamp 
Damped
Sine function (exponential decay) 
1D  Amplitude Phase_Shift Period
BackGround Decay 
strline 
Straight
line 
1D  Gradient Background 
triangl 
Triangular

1D  Amplitude Centre HalfWidth Background 
tophat 
TopHat
rectangular function 
1D  Amplitude Centre HalfWidth Background 
twoexp 
Two exponential decay
functions 
1D  Amplitude1 Tau1 Amplitude2
Tau2 Background 
voigt 
Voigt
function 
1D  Amplitude Centre
HalfWidth_Gauss HalfWidth_Lorz Background 
gauss2d 
Gaussian
function with tilt angle 
2D  Amplitude Centre_X Center_Y
HalfWidth_X HalfWidth_Y Angle Background 
lorz2d 
Lorentzian function with tilt angle  2D  Amplitude Centre_X Center_Y HalfWidth_X HalfWidth_Y Angle Background 
plane2d 
Planar function 
2D  Slope_X Slope_Y Background 
pseudovoigt2d 
Pseudo Voigt with tilt angle  2D  Amplitude Centre_X Center_Y
HalfWidth_X HalfWidth_Y Angle Background LorentzianRatio 
quad2d 
Quadratic (parabola) with tilt angle  2D  Amplitude Centre_X Center_Y Curvature_X Curvature_Y Angle Background 
gaussnd 
ndimensional
Gaussian 
nD 

sf_hard_spheres 
Hard Sphere structure factor
[PercusYevick] 
1D  R rho 
rietveld 
Rietveld
refinement of powder sample with full McStas instrument model 
1D,
2D, 3D 
sample structure, instrument parameters 
sqw_sine3d  Phonon dispersions as sine wave in HKL with a damped harmonic oscillator energy dispersion  4D (HKLw) 
zone center, energy gaps, periodicity, DHO
width, temperature, background 
sqw_spinw 
Spinwave dispersion in HKL using SpinW.  4D (HKLw) 
energy broadening, Temperature, Amplitude,
coupling parameters J... 
sqw_vaks  Phonon dispersions in perovskite cubic crystals using the Vaks parameterisation  4D (HKLw) 
acoustic and optical energies, coupling parameters, soft mode frequency, DHO width, temperature, background 
sqw_cubic_monoatomic  Phonon dispersions in a monoatomic cubic crystal using the Dynamical matrix.  4D (HKLw) 
acoustic force constant ratio and energy scaling, DHO width, temperature, background 
sqw_phonons  Phonon dispersions from the Dynamical matrix, using forces estimated by abinitio DFT using ASE or PHON/QE.  4D (HKLw) 
Creation: POSCAR,
CIF,
or PDB, ... Then, only the DHO line shape. abinitio implies no (few) tunable parameter. 
sqw_linquad 
A phenomenological dispersion which can
describe e.g. an acoustic/linear mode, with quadratic expansion in other directions. This model can be
considered as a local expansion in series of any dispersion. 
4D (HKLw) 
Energy and location of 'gap', slopes, directions of slopes, DHO width, temperature, background 
sqw_acoustopt  A phenomenological dispersion which can describe e.g. a pure acoustic or optical mode, with quadratic expansion around a minimum.  4D (HKLw)  Energy and location of 'gap', slopes, directions of slopes, DHO width, temperature, background 
>> gauss
>> a = gauss;
The list of all available fit functions can be obtained from the command:
>> fits(iData);which also produces the plot above.
We list below a number of models used to describe neutron and
xray scattering from matter. The structure factor S(q) accounts
for the structure of matter at the atomistic/molecular level,
whereas the form factor P(q) accounts for the geometrical
arrangement of large scale scattering units (micelles, tubes,
...). In practice, the scattering from a material can be described
by:
I(q)
= P(q).S(q)
where q is the momentum exchange in the material.
These models have been extracted from:
Form factors: small angle 
Description 
Dimensionality 
Parameters 
ff_core_shell 
Spherical/core shell form
factor [Guinier] 
1D 
R1 R2 eta1 eta2 
ff_sphere

Sphere form factor [Guinier] 
1D 
R eta 
The rietveld model performs a structure refinement (atom
type and position, structure group) of a powder by comparing a
measured diffractogram with a simulated diffractogram using McStas and CrysFML.
This model requires external software to be installed on your
computer. See Requirements below.
Powder diffraction 
Description 
Dimensionality 
Parameters 
rietveld 
Rietveld
refinement of powder sample with full McStas instrument model 
1D,
2D, 3D 
sample structure, instrument parameters 
The 'rietveld' model allows to prepare a sample+instrument model in order to fit a diffraction/structure data set.
>> model=rietveld(structure, ...., instrument, ....)
>> p = fits(model, measurement);
Sample.cell = [10.242696 10.242696 10.242696 90.000 90.000 90.000]; Sample.Spgr = 'I 21 3'; Sample.Ca1 = [0.46737 0.00000 0.25000 0.60046 0.50000 0.0 2.0];
...
>> model = rietveld([ ifitpath 'Data/Na2Ca3Al2F14.cfl' ], 'templateDIFF.instr', 'Powder=reflections.laz; lambda="2.36"; monitors=BananaTheta');Once built, it is possible to set constraints on the model with the syntax such as (see iFunc page) :
>> model.parameter='fix' % to lock its value during a fit processThen we import a data set
>> model.parameter='clear' % to unlock value during a fit process >> model.parameter=[min max] % to bound value
>> model.parameter=[nan nan] % to remove bound constraint
>> model.parameter='' % to remove all constraints on 'parameter'
>> model.Constraint='' % to remove all constraints
>> measurement = iData([ ifitpath 'Data/nac_1645179.dat' ]);The refinement is then obtained by starting :
>> parameters = fits(model, measurement)
with optional arguments as described in the Fit
page (the 'constraints' argument of fits is partly redundant with the model
constraints seen above).The CrysFML does not need to be installed, as the used bits are assembled in the cif2hkl fortran programme, which generates HKL dF2 reflection lists suitable for the PowderN and Isotropic_Sqw McStas components. The cif2hkl programme is part of the iFit distribution, and is compiled by calling the local fortran compiler (gfortran) on demand.cd /etc/apt/sources.list.d sudo wget http://packages.mccode.org/debian/mccode.list sudo aptget update
sudo aptget install mcstassuite
The sqw_sine3d model provides a simple way to model phonontype
dispersions, including simple spinwaves,
acoustic and optical modes, incommensurate dispersions.
Limitation: this model only handles simple sine dispersions, and can not treat mode exchange (interferences). For more advanced spinwave models, use the sqw_spinwave Model (below).
Each dispersion is a sine wave which goes continuously from
energy E0 to E1, along 3 principal lattice directions (HKL). The
dispersion has an energy width (DHO). Schematically, the
dispersion relation is:
w(Q) = E0 + (E1E0)*sin(Q_freq*pi*(QQ0));
along principal axes
where the wavevector/momentum Q is expressed in reciprocal
lattice units [r.l.u]. The parameters of this model allow extended
flexibility in the description of the mode. Along the 3D HKL
volume, a dispersion is described with 10 parameters.
The Q_freq parameter indicates how many dispersion sine
'arches' there are per reciprocal lattice unit [rlu]. A Q_freq
of 1/2 means the dispersion extends from e.g. Q=0 to Q=2 rlu. A Q_freq
of 1 means it extends from Q=0 to Q=1 rlu, and a Q_freq of
2 means there are two arches between Q=0 and Q=1 rlu (all these
with Q0=0). A Q_freq of 0 sets a flat dispersion.
To create the model without defined parameter values, you may
use:
or alternatively, to define starting parameters:sw = sqw_sine3d;
sw = sqw_sine3d(p);where p=[...] is a vector containing the parameter values. It may be given as 1,2,3 and 14 value vector, as detailed below.
A spinwave could for instance mostly use Q0=0,
Q_freq=1, E0=0, E1>0 (2 arches from Q=0 to 1 rlu). In a
simple antiferromagnet, the gap width is E1E0=4J.S with
J=exchange energy and S=magnetic moment of spins.
>> sw = sqw_sine3d([ E0 E1 Q_freq ]) % creates a dispersion from E0 to E1 with given Q frequency, e.g. .5, 1 or 2
A phonon acoustic branch could use Q0=0, Q_freq=.5, E0=0
(1 arch from Q=0 to 1 rlu).
>> acoustic = sqw_sine3d(Emax) % creates an acoustic dispersion up to Emax
An phonon optical branch could use Q0=0, Q_freq=.2, E0>E1 E1>0 (1 arch from Q=0 to 1 rlu with Q=0 energy  Raman frequency).
>> optical = sqw_sine3d([ E0 E1 ]) % creates an optical dispersion from E0 to E1
The model parameters allow to tune the dispersion:
The model parameters are:
p( 1)= E1_qh energy at QH half period [meV]
p( 2)= E1_qk energy at QK half period [meV]
p( 3)= E1_ql energy at QL half period [meV]
p( 4)= E0 zonecentre energy gap [meV]
p( 5)= QH0 QH zonecentre [rlu]
p( 6)= QK0 QK zonecentre [rlu]
p( 7)= QL0 QL zonecentre [rlu]
p( 8)= QH_freq QH frequency [multiples of pi]
p( 9)= QK_freq QK frequency [multiples of pi]
p( 10)= QL_freq QL frequency [multiples of pi]
p( 11)= Gamma Damped Harmonic Oscillator width in energy [meV]
p( 12)= Temperature [K]
p( 13)= Amplitude
p( 14)= Background
The axes needed for the evaluation are expressed in rlu for
QH,QK,QL and in meV for the energy.
The Model evaluation is as usual (in 4D):
model(p, qh, qk, ql, w) % return a matrix. p can be [] to use guessed/default parameters
iData(model, p, qh, qk, ql, w) % return an iData.
A usage example is as follows:
>> ac=sqw_sine3d(5); % an acoustic branch up to 5 meV
>> qh=linspace(0,1,50);qk=qh; ql=qh'; w=linspace(0.01,10,50); % the axes for evaluation
>> f=iData(ac,[],qh,qk,ql,w); % evaluate the model onto given axes
>> plot3(log(f(:,:,1,:))); % plot as volume in [QH,QK,w, QL=0]. You can also use surf and scatter3 for other rendering
The axes are given as vectors, but the second is made
nonparallel to the others, to indicate we wish to build a volume
out of this. Without transposing the vector, as all axes are the
same length, they would be interpreted as event data, and the
resulting evaluation would only contain 50 values.
It is possible to stack as many modes as possible, in different
flavors. In this case it is advisable to link e.g. the Temperature
and Background parameters:
>> acoustic = sqw_sine3d(5); optical = sqw_sine3d([ 10 8 ]); sw = sqw_sine3d([ 2 4 1 ]);
>> disp3 = acoustic + optical + sw;
>> disp3.Temperature_2 = '"Temperature"'; % Temperature_2 = Temperature (1st sin3d)
>> disp3.Temperature_3 = '"Temperature"'; % Temperature_3 = Temperature
>> disp3.Name='sqw_sine3d: acousitc+optical+sw';
>> mlock(disp3, {'Background_2','Background_3'}); % keep them as 0 (default)
>> qh=linspace(0,1,50);qk=qh; ql=qh; w=linspace(0.01,10,51); % the axes for evaluation
>> f=iData(disp3,[],qh,qk,ql,w); % evaluate using initial model parameters
>> plot3(log(f(:,:,1,:))); % plot as volume in [QH,QK,w, QL=0].
The plot (surf), plot3, scatter3, and slice
methods for plotting all provide nice looking rendering of volume
data. See Plot/3D page.
The sqw_spinw model can embed a spinwave model from the
SpinW to "
simulate magnetic structures and excitations of given spin
Hamiltonian using classical Monte Carlo simulation and linear
spin wave theory."
Basically, SpinW allows to model a crystal structure including a
spin Hamiltonian in the form:
S(q,w) 
Description 
Dimensionality 
Parameters 
sqw_spinw 
Spinwave dispersion in HKL
using SpinW. 
4D
(HKLw) 
energy broadening,
Temperature, Amplitude, parameters J... Axes are in rlu. 
s = sqw_spinw(sw);or
s = iFunc(sw);
The SpinW toolbox allows
to build spinwave models incuding spinspin pair coupling,
anisotropy terms and external field contribution. The model is
then converted into an iFunc object, so that fitting and plotting
is straightforward. An easy way to build a SpinW object is to read
a CIF
file, with syntax:
s = sqw_spinw(sw('cif_file'))
Additional options can be used when creating the Model, with:
where 'options' is a structure that can hold:s = sqw_spinw(sw, options);
p(1)= Gamma energy broadening (instrumental) [meV]The axes needed for the evaluation are expressed in rlu for QH,QK,QL and in meV for the energy [1 meV = 241.8 GHz = 11.604 K = 0.0965 kJ/mol]. The model can be evaluated by giving axes and model parameters (see examples below).
p(2)= Temperature of the material [K]
p(3)= Amplitude
p(4...)= coupling parameters of the Hamiltonian (J)
model(p, qh, qk, ql, w) % return a matrix. p can be [] to use guessed/default parametersAs an example, we plot a simple square lattice Heisenberg antiferromagnet model:
iData(model, p, qh, qk, ql, w) % return an iData
sq = sw_model('squareAF',2,0); % create the SW object (using SpinW)References
s=sqw_spinw(sq); % create the Model
qh=linspace(0.01,1.5,50);qk=qh; ql=qh'; w=linspace(0.01,10,50);
f=iData(s,s.p,qh,qk,ql,w); plot3(log(f(:,:,1,:))); % evaluate and plot
S. Toth and B. Lake, J. Phys.: Condens. Matter 27, 166002 (2015).Requirements
In principle, the SpinW Matlab toolbox must have been installed from e.g. <https://github.com/tsdev/spinw> or <https://www.psi.ch/spinw/spinw>.
The version 2.1 of SpinW is included in iFit, so that the sqw_spinw Model does not require any additional installation.
The sqw_vaks model computes the dispersion of the 5
lowest phonon
dispersions in perovskite cubic crystals. It is based on the Vaks
parametrization (see references below). The TA1,TA2,LA,TO1 and TO2
dispersions are obtained from 8 parameters, and a DHO line shape
is added. The dispersions are anisotropic, with 'valley' and soft
mode.
The dynamical matrix is 5x5 and its eigenvalues are the mode
frequencies.
Limitation: even though this is model uses
few parameters, the dynamic range is limited to e.g. q
< 0.30.5 rlu and w < 100 meV.
To create the model without defined parameter values, you may
use:
or alternatively, to define starting parameters:s = sqw_vaks;
s = sqw_vaks(p);where p=[...] is a vector containing the parameter values. It may be given as a 12 value vector or a string, as detailed below.
p( 1)= At transverse acoustic slope [meV^{2}/rlu^{2}]The axes needed for the evaluation are expressed in rlu for QH,QK,QL and in meV for the energy [1 meV = 241.8 GHz = 11.604 K = 0.0965 kJ/mol]. The model can be evaluated by giving axes and model parameters (see examples below).
p( 2)= Al longitudinal acoustic slope [meV^{2}/rlu^{2}]
p( 3)= Aa anisotropic acoustic slope [meV^{2}/rlu^{2}]
p( 4)= St soft mode transverse soft optical slope [meV^{2}/rlu^{2}]
p( 5)= Sa soft mode anisotropic soft optical slope [meV^{2}/rlu^{2}]
p( 6)= Vt transverse acousticoptical coupling [meV^{2}/rlu^{2}]
p( 7)= Va entered as: anisotropic acousticoptical coupling [meV^{2}/rlu^{2}]
p( 8)= w0 soft mode frequency at q=0, depends on temperature [meV]
p( 9)= Gamma Damped Harmonic Oscillator width in energy [meV]
p( 10)= Temperature [K]
p( 11)= Amplitude
p( 12)= Background
>> s=sqw_vaks('KTaO3'); % create model, with KTaO3 parameters
>> qh=linspace(0,.5,50);qk=qh; ql=qh'; w=linspace(0.01,10,51);
>> f=iData(s,[],qh,qk,ql,w); % evaluate model into an iData
>> scatter3(log(f(:,:,1,:)),'filled'); % plot ql(1)=0 plane
w(k) = E_{max}*sin(k.a) along principal axes.The dynamical matrix is 3x3 with sine and cosine terms. Eigenvalues provide the mode energies.
p( 1)= C_ratio C1/C2 force constant ratio first/second neighboursThe axes needed for the evaluation are expressed in rlu for QH,QK,QL and in meV for the energy [1 meV = 241.8 GHz = 11.604 K = 0.0965 kJ/mol]. The model can be evaluated by giving axes and model parameters (see examples below).
p( 2)= E0 sqrt(C1/m) energy [meV]
p( 3)= Gamma Dampled Harmonic Oscillator width in energy [meV]
p( 4)= Temperature [K]
p( 5)= Amplitude
p( 6)= Background
>> s=sqw_cubic_monoatomic([ 3 3 ]); % create model, with C1/C2=3 E0=3The plot (surf), plot3, scatter3, and slice methods for plotting all provide nice looking rendering of volume data. See Plot/3D page.
>> qh=linspace(0,.5,50);qk=qh; ql=qh'; w=linspace(0.01,10,51);
>> f=iData(s,[],qh,qk,ql,w); % evaluate model into an iData
>> plot3(log(f(:,:,1,:))); % plot ql(1)=0 plane
S(q,w) 
Description 
Dimensionality 
Parameters 
sqw_linquad 
A phenomenological dispersion which can describe an acoustic or optical mode. This model can be considered as a local expansion in series of any dispersion.  4D (HKLw) 
Energy and location of 'gap', slopes, directions of slopes, DHO width, temperature, background 
p( 1)= DC_Hdir1 Slope1 dispersion direction, linear, H [rlu]
p( 2)= DC_Kdir1 Slope1 dispersion direction, K [rlu]
p( 3)= DC_Ldir1 Slope1 dispersion direction, L [rlu]
p( 4)= DC_Hdir2 Slope2 dispersion direction, H (transverse) [rlu]
p( 5)= DC_Kdir2 Slope2 dispersion direction, K (transverse) [rlu]
p( 6)= DC_Ldir2 Slope2 dispersion direction, L (transverse) [rlu]
p( 7)= DC_Slope1 Dispersion slope along 1st axis, linear [meV/rlu]
p( 8)= DC_Slope2 Dispersion slope along 2nd axis (transverse to 1st, in plane) [meV/rlu]
p( 9)= DC_Slope3 Dispersion slope along 3rd axis (transverse to 1st, vertical) [meV/rlu]
p( 10)= Ex_H0 Excitation location H [rlu]
p( 11)= Ex_K0 Excitation location K [rlu]
p( 12)= Ex_L0 Excitation location L [rlu]
p( 13)= Ex_E0_Center Excitation location, Energy [meV]
p( 14)= DHO_Amplitude
p( 15)= DHO_Damping Excitation damping, halfwidth [meV]
p( 16)= DHO_Temperature Temperature [K]
p( 17)= Background
velocity_{[m/s]} = Slope_{[meV/rlu]} *a_{[Angs]} /2pi * 151.9A usage example is:
This model is suited to describe a quadratic dispersion with given minimum/gap.
The dispersion slopes (in [meV/rlu]) are given for user defined axes
(e.g. principal crystal directions in reciprocal space). Two
directions must be given, which are used to define an orthonormal
reciprocal coordinate frame. The dispersion slopes apply
into this frame.
The general dispersion relation has the form:
S(q,w) 
Description 
Dimensionality 
Parameters 
sqw_acoustopt 
A phenomenological dispersion which can
describe an acoustic or optical mode. 
4D (HKLw) 
Energy and location of 'gap', slopes, directions of slopes, DHO width, temperature, background 
p( 1)= DC_Hdir1 Slope1 dispersion direction, H [rlu]The Slope parameters allow to estimate the mode phase velocity for acoustic modes when E0=0. Using the definition rlu = 2pi/a where a is the excitation periodicity, we can derive:
p( 2)= DC_Kdir1 Slope1 dispersion direction, K [rlu]
p( 3)= DC_Ldir1 Slope1 dispersion direction, L [rlu]
p( 4)= DC_Hdir2 Slope2 dispersion direction, H (transverse) [rlu]
p( 5)= DC_Kdir2 Slope2 dispersion direction, K (transverse) [rlu]
p( 6)= DC_Ldir2 Slope2 dispersion direction, L (transverse) [rlu]
p( 7)= DC_Slope1 Dispersion slope along 1st axis [meV/rlu]
p( 8)= DC_Slope2 Dispersion slope along 2nd axis (transverse to 1st, in plane) [meV/rlu]
p( 9)= DC_Slope3 Dispersion slope along 3rd axis (transverse to 1st, vertical) [meV/rlu]
p( 10)= Ex_H0 Minimum of the dispersion, H [rlu]
p( 11)= Ex_K0 Minimum of the dispersion, K [rlu]
p( 12)= Ex_L0 Minimum of the dispersion, L [rlu]
p( 13)= Ex_E0_Center Minimum of the dispersion, Energy [meV]
p( 14)= DHO_Amplitude
p( 15)= DHO_Damping Excitation damping, halfwidth [meV]
p( 16)= DHO_Temperature Temperature [K]
p( 17)= Background
velocity_{[m/s]} = Slope_{[meV/rlu]} *a_{[Angs]} /2pi * 151.9A usage example is:
The sqw_phonons model computes phonon dispersions using
the Atomic Simulation
Environment (ASE) in the socalled small displacement
methodology. The only needed input for the model to run is the
name of a configuration file describing the material lattice and
atom positions. This file can be in any ASE
supported format, e.g. POSCAR,
CIF,
PDB, ... This model
requires to install external software on the computer. See Requirements below.
S(q,w) 
Description 
Dimensionality 
Parameters 
sqw_phonons 
Phonon dispersions from the Dynamical matrix, using forces estimated by abinitio using ASE and a selection of DFT codes (EMT, GPAW, ABINIT, Elk, Dacapo, NWChem, QuantumEspresso, VASP).  4D (HKLw) 
Creation: POSCAR,
CIF,
PDB, ... Then, only the DHO line shape. abinitio implies no (few) tunable parameter. Axes in rlu. 
D_{α}_{β}=1/√(M_{α}M_{β}) ∂^{2}V/(∂u_{α}∂u_{β}) where V is the interatomic potential, M is the atom mass, αβ are atom indices, u is the displacement of these atoms around their equilibrium position._{ }The equation of motions is then:
d^{2}u_{α}/dt^{2} =  ∑ D_{α}_{β}.u_{β}In order to simplify this expression, we assume a periodic motion:
u_{α}=U_{α}(k) e ^{i(k.r  ωt}^{)}where r is the position of the atom in a socalled supercell, which replicates the lattice primitive cell, and its extent sets the potential extension distance, which is assumed to vanish at long distance. k and ω are the momentum and energy of the wave.
s = sqw_phonons; % displays an entry dialog box and proceeds with a fully automatic computation, and plots final results. With QuantumEspresso, a dialog may show up to select pseudopotentials when there is more than one possibility.The general syntax for creating the S(q,w) model is any of:
s = sqw_phonons('defaults'); % creates a simple fcc Al model and QE as calculator
s = sqw_phonons(configuration);
s = sqw_phonons(configuration, options, ...);The result is a 4D S(hkl,w) phonon dispersion Model (iFunc object).
s = sqw_phonons(configuration, calculator, options, ...);
>> save 'my_model' modelYou can later load this model again within iFit, without the need to recompute forces (which may be a long process with DFT codes).
>> load 'my_model'
The model parameters, once built, are the following:
p( 1)= AmplitudeThe axes needed for the evaluation are expressed in rlu for QH,QK,QL and in meV for the energy [1 meV = 241.8 GHz = 11.604 K = 0.0965 kJ/mol]. The model can be evaluated by giving axes and model parameters (see examples below).
p( 2)= Gamma Damped Harmonic Oscillator width in energy [meV]
p( 3)= Background
p( 4)= Temperature [K]
>> qh=linspace(0.01,.5,50);qk=qh; ql=qh; w=linspace(0.01,50,51);
>> f=iData(s,[],qh,qk,ql',w); % evaluate model into an iData.
>> plot3(log(f(1,:,:,:))); % plot QH=0 plane
configuration 
Description 
a file name [ ifitpath 'Data/POSCAR_Al'] 
Any file in a format supported by ASE,
including CIF, ETSF, Gromacs, LAMMPS, PDB, SHELX,
VASP/POSCAR, MaterialStudio, XYZ. See full
list. You can use as example (provided in iFit Data
directory):POSCAR_Al 
'bulk("Cu",
"fcc", a=3.6, cubic=True)' 
A bulk material, with formula (Cu, MgO,
NaCl), lattice constant and structure (sc, fcc, bcc, hcp,
diamond, zincblende, rocksalt, cesiumchloride, fluorite or
wurtzite). See description of bulk.bulk("Cu", "fcc", a=3.6) 
'molecule("isobutane")' 
Any molecule
known to ASE, e.g. 'isobutene',
'CH3CH2OH', 'CH3COOH', 'COF2', 'CH3NO2', 'CF3CN',
'CH3OH', 'CCH', 'CH3CH2NH2', 'PH3', 'Si2H6', 'O3', 'O2',
'BCl3', 'CH2_s1A1d', 'Be', 'H2CCl2', 'C3H9C', 'C3H9N',
'CH3CH2OCH3', 'BF3', 'CH3', 'CH4', 'S2', 'C2H6CHOH',
'SiH2_s1A1d', 'H3CNH2', 'CH3O', 'H', 'BeH', 'P',
'C3H4_C3v', 'C2F4', 'OH', 'methylenecyclopropane',
'F2O', 'SiCl4', 'HCF3', 'HCCl3', 'C3H7', 'CH3CH2O',
'AlF3', 'CH2NHCH2', 'SiH2_s3B1d', 'H2CF2', 'SiF4',
'H2CCO', 'PH2', 'OCS', 'HF', 'NO2', 'SH2', 'C3H4_C2v',
'H2O2', 'CH3CH2Cl', 'isobutane', 'CH3COF', 'HCOOH',
'CH3ONO', 'C5H8', '2butyne', 'SH', 'NF3', 'HOCl',
'CS2', 'P2', 'C', 'CH3S', 'O', 'C4H4S', 'S', 'C3H7Cl',
'H2CCHCl', 'C2H6', 'CH3CHO', 'C2H4', 'HCN', 'C2H2',
'C2Cl4', 'bicyclobutane', 'H2', 'C6H6', 'N2H4',
'C4H4NH', 'H2CCHCN', 'H2CCHF', 'cyclobutane', 'HCl',
'CH3OCH3', 'Li2', 'Na', 'CH3SiH3', 'NaCl', 'CH3CH2SH',
'OCHCHO', 'SiH4', 'C2H5', 'SiH3', 'NH', 'ClO', 'AlCl3',
'CCl4', 'NO', 'C2H3', 'ClF', 'HCO', 'CH3CONH2',
'CH2SCH2', 'CH3COCH3', 'C3H4_D2d', 'CH', 'CO', 'CN',
'F', 'CH3COCl', 'N', 'CH3Cl', 'Si', 'C3H8', 'CS', 'N2',
'Cl2', 'NCCN', 'F2', 'CO2', 'Cl', 'CH2OCH2', 'H2O',
'CH3CO', 'SO', 'HCOOCH3', 'butadiene', 'ClF3', 'Li',
'PF3', 'B', 'CH3SH', 'CF4', 'C3H6_Cs', 'C2H6NH', 'N2O',
'LiF', 'H2COH', 'cyclobutene', 'LiH', 'SiO', 'Si2',
'C2H6SO', 'C5H5N', 'transbutane', 'Na2', 'C4H4O',
'SO2', 'NH3', 'NH2', 'CH2_s3B1d', 'ClNO', 'C3H6_D3h',
'Al', 'CH3SCH3', 'H2CO', 'CH3CN'molecule("H2O") 
'nanotube(6,
0, length=4)' 
A carbon (or other material) nanotube. See
description of nanotube. 
'crystal(...)' 
The generic crystal builder from ASE. See
description here. NaCl: Rutile:crystal(["Na", "Cl"], [(0, 0, 0), (0.5, 0.5, 0.5)], spacegroup=225, cellpar=[5.64,5.64,5.64, 90, 90, 90]) crystal(["Ti", "O"], basis=[(0, 0, 0), (0.3, 0.3, 0.0)], spacegroup=136, cellpar=[4.6,4.6,2.95, 90, 90, 90]) 
>> s=sqw_phonons(POSCAR,'ABINIT','nsteps=100; occupations=smearing; ecut=1000; pps=paw;');More options can be obtained from 'help sqw_phonons'
The pseudopotentials are given as many
filenames/directories (UPF format) suitable for QuantumEspresso.
When not given, available potentials will be searched in current
location, the Home directory, and the standard QuantumEspresso
locations (PSEUDO_DIR, $HOME/espresso/pseudo). When more than one
potential is available for an atom type, a list selector is
displayed to request which one to choose. The PBEPAW potentials
are preferred and preselected. Pseudopotential can be obtained
at QuantomEspresso
PP selector. We recommend the SSSP library, which is
very accurate.
If you have difficulties in the model SCF convergence, you may use e.g:
>> s=sqw_phonons(POSCAR,'QuantumEspresso', 'random','nsteps=200; toldfe=1e6; occupations=smearing; ecut=500; mixing_ndim=12; miximg_beta=0.3;');If you are tight on the memory (e.g. the system is large), you may try (but it may not converge):
>> s=sqw_phonons(POSCAR,'QuantumEspresso', 'random','ecut=15; mixing_ndim=4; diagonalization=cg');
The sqw_phonons model may use any of the calculators
listed below:
'calculator=...' 
Description 
time elapsed [s] 
Delta [meV/atom] 
EMT  Effective Medium Theory calculator (only
for Al,Cu,Ag,Au,Ni,Pd,Pt,H,C,N,O). Very fast, but
limited. 
0.9 

QuantumEspresso  Densityfunctional theory, plane waves, and
pseudopotentials based code for electronicstructure
calculations and materials modeling. This options uses PHON and
Quantum Espresso.
Our best choice with the SSSP
library pseudo. This calculator can not compute the mode intensity/polarisation vectors. 
319 
SSSP: 0.3 
QuantumEspresso_ASE  Same as above, but using the ASE interface
instead of PHON.
The computation is about 36 times slower as when using PHON, as
it often uses more atom displacements (with accurate setting). 
957 
SSSP: 0.3 
ABINIT  Planewave pseudopotential code. Can be used
with JTH
and GBRV
pseudo lib (nsteps=100 ecut=1000 eV) 
1742 (7.10.5) 
pps=JTH: 0.4 pps=GBRV: 0.8 pps=HGH.k: 2.1 
VASP 
Planewave PAW code. Requires a
specific license. 
3741 (4.6) 

GPAW 
Realspace/planewave/LCAO PAW code.  4548 
1.5 
ELK  Full PotentialAll Electrons LAPW code. A
very efficient allelectrons code. 
4788 
0.3 
Dacapo  Planewave ultrasoft pseudopotential code. Not recommended for phonon calculations.  
NWChem  Gaussian based electronic structure code. Not
recommended for phonon calculations. 
It is generally recommended to use the following settings for
very accurate phonon band computations:
If the SCF convergence fails, you may try:
>> s=sqw_phonons(POSCAR, 'nsteps=200; toldfe=1e6; occupations=0.27; ecut=1500; ');For insulators, you may also use a very small smearing, e.g. 0.01 eV.
>> s=sqw_phonons(POSCAR, 'ecut=500; diagonalization=cg');
>> s=sqw_phonons([ ifitpath 'Data/POSCAR_Al'],'metal','EMT');The same procedure can be applied on an evaluated iData object out of the model.
>> pow=sqw_powder(s); plot(log(pow));
>> qh=linspace(0.01,3,50);qk=qh; ql=qh; w=linspace(0.01,50,51);In both cases, when the reciprocal lattice information is found in the object (matrix B=[a* b* c*]), the conversion from 'rlu' (e.g. 2pi/a units) to 'Angs^{1}' is automatically done (and this information is stored in the Model upon creation).
>> f=iData(s,[],qh,qk,ql,w);
>> pow=sqw_powder(f); plot(log(pow));
>> s=sqw_phonons('bulk("Cu", "fcc", a=3.6)','EMT','metal'); % create model from bulkTo get a system with acoustic and optical modes:
>> qh=linspace(0.01,.5,50);qk=qh; ql=qh; w=linspace(0.01,50,51);
>> f=iData(s,[],qh,qk,ql,w); % evaluate model into an iData
>> plot3(log(f(:,:,1,:))); % plot ql=0 plane
>> s=sqw_phonons([ ifitpath 'Data/POSCAR_AlN'],'EMT'); % create model from given POSCARIn addition, the phonon Density of States can be computed, when specifying the 'dos' options upon model creation. This computation can be rather intensive. It is then computed during the first model evaluation (subsequent evaluations will be faster), and stored in the model UserData.DOS so that it can be viewed with e.g. :
>> qh=linspace(0.01,.5,50);qk=qh; ql=qh; w=linspace(0.01,50,51);
>> f=iData(s,[],qh,qk,ql,w); % evaluate model into an iData
>> scatter3(log(f(:,:,1,:))); % plot ql=0 plane with acoustic and optical modes
>> s=sqw_phonons([ ifitpath 'Data/POSCAR_Al'],'dos=1','EMT'); % create model from given POSCAR, with DOS computationThe plot (surf), plot3, scatter3, and slice methods for plotting all provide nice looking rendering of volume data. See Plot/3D page.
>> qh=linspace(0.01,.5,50);qk=qh; ql=qh; w=linspace(0.01,50,51);
>> f=iData(s,[],qh,qk,ql,w); % evaluate model into an iData
>> figure; plot(s.UserData.DOS);
External calculators should be installed separately. For Debian systems, we provide most packages at packages.mccode.org. Type (Ubuntu/Mint/Debian):% sudo aptget install pythonase
% sudo aptaddrepository 'deb http://packages.mccode.org/debian stable main'
% sudo aptget update
% sudo aptget install ifit elklapw dacapo dacapopsp gpaw gpawsetups nwchem nwchemdata pythonmatplotlib
% sudo aptget install abinit abinitdoc abinitpsp quantumespresso quantumespressosssp phon libnetcdfdev libxc1 openmpibin
% sudo aptget install abinitpspjth abinitpspgbrv pythonaseqespglib
Note: The gpawsetups package is named gpawdata in Ubuntu distributions >= 15.10.Pseudo potentials:
>> h=ifitmakefunc;The only required argument is the expression, but you may as well enter a function name and description, parameter names and default values, as well as a constraint to be evaluated before computing the function value. When started without any parameters, the Gaussian function settings are used.
>> h=edit(iFunc); % equivalent
>> [p,criteria,message,output]= fits(a, h);A quick definition, without using the dialog, can be performed with the syntax:
>> [p,criteria,message,output]= fits(h, a); % same as above
>> h=ifitmakefunc('p(1)*exp( (xp(2))/p(3) )');where p is the vector that holds parameters. Axes are x (rows),y (columns),z (pages),t (time). Beware to use a model expression with elementwise division and multiplication operators (./ and .*). The function builder is known to work well for 1D and 2D functions.
>> fits(a, iFunc('p(1)*exp( (xp(2))/p(3) )')); % same as above
>> h = gauss + lorzThe full syntax for the model builder is:
>> h = convn(lorz, 3)
>> h = convn(gauss, lorz)
>> h = gauss .* lorz; % a 1D model
>> h = gauss * lorz; % a 2D model
>> h=ifitmakefunc('expression');which is equivalent to:
>> h=ifitmakefunc(iFunc model);
>> h=edit(iFunc model); % same as above
>> h=ifitmakefunc('function_name', 'description', 'Parameter1 Parameter2 ...', ...
'expression', [default parameter values], 'constraint');
>> fun.Name='function short name';The generated code will basically be:
>> fun.Description='function description'; % function long description
>> fun.Parameters='Parameter1 Parameter2 ...'; % parameter names. When empty, names are given according to the Expression analysis (when appropriate)
>> fun.Guess=[0 1 2 ...] or 'automatic'; % parameter default values (vector), or automatic mode
>> fun.Expression='expression using p and x,y,z,t...'; % value of the model (required). p is the vector that holds parameters. Axes are x,y,z,t,u,v,w.
>> fun.Constraint='evaluated before Expression'; % constraint evaluated prior to the model Expression
>> h=iFunc(fun);
>> fits(a, h);
function y=Name(p, axes, ...)The resulting function has the ability to identify itself (disp(model) provides detailed informations), compute automatic starting parameters with e.g. feval(model, 'guess'), display itself plot(model), and evaluate its value of course feval(model, parameter_values, axes, ...). The resulting model is a iFunc object. Refer to that object description for more information.
% Description
% Parameters: [fun.Parameters]
Constraint;
y=Expression;
>> h=gauss+lorz; h.Constraint = 'p(8)=0';which creates a new function which is the sum of a Gaussian and a Lorentzian. Unspecified arguments are guessed/automatically set. The second redundant Lorentzian Background p(8) parameter is forced to 0 so that it does not correlate with the Gaussian Background p(4).
>> c = fconv(a,b); % returned convoluted object with size which is size(a)+size(b)+1When some data has to be convoluted with a response function (e.g. instrument resolution function), the usual options to use should be:
>> c = fconv(a,b, 'same'); % returned convoluted object with size which is size(a)
>> c = fconv(a,b, 'valid'); % returned convoluted object with size which is size(a)size(b)+1
>> c = fconv(a,b, 'pad'); % pads 'a' with starting/ending values to minimize border effects
>> c = fconv(a,b, 'center'); % centers 'b' so that convolution does not shift 'a'
>> c = fconv(a,b, 'norm'); % normalizes 'b' so that convolution does not change 'a' integral
>> c = fconv(a,b, 'background');% subtracts minimal value in 'b' so that convolution does not change 'a' integral
>> c = fconv(a,b, 'deconv'); % deconvolution, but very sensitive to noise (use with caution)
>> c = fconv(a,b, 'same pad background center norm');These convolution methods have been ported to iFunc models as
>> c = fconvn(a,b); % same as above in a shorter call
>> a = convn(lorz, 3) % convolution of a Lorentzian with a Gaussian of width 3where the vector/matrix 'b' holds the response function (filter) with the same axis binning as the object 'a'.
>> a = convn(lorz, gauss) % a Voigt function...
>> a = convn(lorz, 'double(b)'); a.Constraint = 'global b'; % convolute with a global variable 'b'
>> h = convn(dho, gauss);
>> h.Constraint='x=linspace(min(x),max(x), 5*length(x))'; % artificially extend the 'x' axis
>> h = h + 'signal=signal(1:5:end);'; % shrink 'signal' back to initial size
>> edit gauss
>> edit voigt
>> ifitmakefunc(voigt) % popup a dialogue
>> edit(voigt) % same as ifitmakefunc for iFunc objects