iFit: Models (fit models)

• How to use models
• List of default fit models
• Specialized models
  • Small Angle Scattering models (neutrons, x-rays)
  • Powder structure refinement (Rietveld)
  • Phonon/spin-wave dispersion (simple)
  • Spin-wave model (SpinW)
  • Spin-wave model (SPINWAVE)
  • Phonon dispersion in perovskites
  • Acoustic phonon dispersion in mono-atomic cubic crystal
  • Phenomenological linear/quadratic dispersion
  • Phenomenological acoustic/optic dispersion (quadratic)
  • Phonon dispersion from ab-initio force estimate and dynamical matrix
• Model builder: defining a new model easy
  • Creating complex functions from simple functions
  • Using signal convolution/correlation in new functions
• How to write manually a model function

How to use models

>> model		% display model information
>> 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;
>> model = gauss+lorz;	% add models
>> 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

	>> iFunc(iData('filename'))

List of default fit models

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	Dose-response curve with variable Hill slope. This is a sigmoid or S-shaped.	1D	Amplitude Center Slope BackGround
expon	Exponential decay	1D	Amplitude Tau Background
expstretched	Exponential - Stretched	1D	Amplitude Tau Exponent Background
gauss	Gaussian where the HalfWidth is in fact σ. The 'true' half width is thus 1.177*HalfWidth.	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	Log-Normal distribution	1D	Amplitude Center Width BackGround
lorz	Lorentzian (aka Cauchy) used with Amplitude uncorrelated to Width.	1D	Amplitude Centre HalfWidth Background
ngauss	multiple Gaussian where the HalfWidth is in fact σ. The 'true' half width is thus 1.177*HalfWidth.	n*1D
nlorz	multiple Lorentzian used with Amplitude uncorrelated to Width.	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 S-shaped 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	Top-Hat 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 where the HalfWidth is in fact σ. The 'true' half width is thus 1.177*HalfWidth.	2D	Amplitude Centre_X Center_Y HalfWidth_X HalfWidth_Y Angle Background
lorz2d	Lorentzian function with tilt angle used with Amplitude uncorrelated to Width.	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	n-dimensional Gaussian	nD
sf_hard_spheres	Hard Sphere structure factor [Percus-Yevick]	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	Spin-wave 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 ab-initio DFT using ASE	4D (HKLw)	Creation: POSCAR, CIF, or PDB, ... Then, only the DHO line shape. ab-initio 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
sqw_spinwave	Spin-wave dispersion in HKL using SPINWAVE (LLB).	4D (HKLw)	Amplitude, and user defined parameters (coupling...)

>> gauss
>> a = gauss;

>> fits(iData);

Specialized models

Small Angle Scattering models (neutrons, x-rays)

We list below a number of models used to describe neutron and x-ray 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:

• SASfit, by Joachim Kohlbrecher and Ingo Bressler PSI (2006-2011)
• Ingo Bressler, Joachim Kohlbrecher and Andreas F. Thunemann, Journal of Applied Crystallography, 2015, 48 (5), 1587-1598

Structure factors	Description	Dimensionality	Parameters
sf_hard_spheres	Hard Sphere structure factor [Percus-Yevick]	1D	R rho
sf_square_well	structure factor of particles interacting with a square well potential [Sharma]	1D	R rho epsilon Delta T
sf_sticky_hard_spheres	Sticky Hard Sphere structure factor [Baxter/Menon]	1D	R rho tau

Powder structure refinement (Rietveld)

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

• CrysFML by Juan Rodriguez-Carvajal and Javier Gonzalez-Platas, ILL and ULL, Tenerife, Spain (LGPL 3.0)
• Commission on Crystallographic Computing, IUCr Newsletter No.1, pp 50-58,  January 2003 [link]
  

>> 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');

>> model.parameter='fix'     % to lock its value during a fit process
>> 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' ]);

>> parameters = fits(model, measurement)

• MaterialsProject<https://www.materialsproject.org/> using VASP. Can export structures as CIF.
  
• Inorganic Crystal Structure Database <http://icsd.ill.fr> with plenty of CIFs.

1. CrysFML, Commission on Crystallographic Computing, IUCr Newsletter No.1, pp 50-58,  January 2003 [link].
2. McStas: K. Lefmann and K. Nielsen, Neutron News 10, 20, (1999) ; P. Willendrup, E. Farhi and K. Lefmann, Physica B, 350 (2004) 735.

cd /etc/apt/sources.list.d
sudo wget http://packages.mccode.org/debian/mccode.list
sudo apt-get update
sudo apt-get install mcstas-suite

Phonon/spin-wave dispersion (simple)

The sqw_sine3d model provides a simple way to model phonon-type dispersions, including simple spin-waves, acoustic and optical modes, incommensurate dispersions.

Limitation: this model only handles simple sine dispersions, and can not treat mode exchange (interferences). For more advanced spin-wave models, use the sqw_spinwave Model (below).

S(q,w)	Description	Dimensionality	Parameters
sqw_sine3d	Phonon dispersions as sine wave in HKL (3D) with a damped harmonic oscillator energy dispersion	4D (HKLw)	zone center, energy gaps, periodicity. Axes are in rlu.

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 + (E1-E0)*sin(Q_freq*pi*(Q-Q0));         along principal axes

where the wave-vector/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:

sw = sqw_sine3d;

sw = sqw_sine3d(p);

A spin-wave 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 anti-ferromagnet, the gap width is E1-E0=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:

• To shift the minimum/maximum Q of the dispersion, move QH0,QK0,QL0 parameters.
• To change the extent of the dispersion in Q, vary QH_freq,QK_freq,QL_freq parameters.
• To change to minimum and maximum energy, move E0 and E1_qh,E1_qk,E1_ql parameters.
• To model an incommensurate dispersion, move both QH0,QK0,QL0 and the QH_freq,QL_freq,QK_freq parameters to incommensurate (non rational) values.
• An anisotropic model can be obtained by having different E1 and Q_freq parameters along axes.
• A 1D or 2D dispersion can be obtained when setting the frequencies (e.g. QK_freq QL_freq) to 0
  

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 zone-centre energy gap [meV]
  p(  5)=                 QH0 QH zone-centre [rlu]
  p(  6)=                 QK0 QK zone-centre [rlu]
  p(  7)=                 QL0 QL zone-centre [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 non-parallel 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.

Spin-wave model (SpinW)

The sqw_spinw model can embed a spin-wave 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	Spin-wave dispersion in HKL using SpinW.	4D (HKLw)	energy broadening, Temperature, Amplitude,  parameters J... Axes are in rlu.

s = sqw_spinw(sw);

s = iFunc(sw);

The SpinW toolbox allows to build spin-wave models incuding spin-spin 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:

s = sqw_spinw(sw, options);

• options.component: a string to specify the component to use as intensity, as documented in sw_egrid function. Default is 'Sperp'. Suggested is also 'Sxx+Syy+Szz'.

p(1)=            Gamma       energy broadening (instrumental) [meV]
p(2)=            Temperature of the material [K]
p(3)=            Amplitude
p(4...)=         coupling parameters of the Hamiltonian (J)

• When all axes are given as vectors of same orientation and same length, the HKLE set is assumed to be list of 'events'.
• When the H K L axes are vectors of same length/orientation, the HKL locations is assumed to be a q-path in the BZ, but the E (energy) axis can be of different length/orientation to compute the Model along (q,E) path, resulting in a 2D object. Typically, the band_structure function does this, to provide the phonon dispersion curve along principal directions.
• When the H K L E axes are given as mixed orientations/length, as vectors or matrices, the Model evaluation is done on a HKLE cube, which is well suited to visualize the whole dispersion surfaces. In the example below, ql has a different orientation (it is transposed) as qh and qk, which triggers a 4D HKLE cube evaluation.

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

sq = sw_model('squareAF',2,0);                  % create the SW object (using SpinW)
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).

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. Updating SpinW may affect the iFit-SpinW interface.

Spin-wave model (SPINWAVE)

model = sqw_spinwave('file');

lmo = sqw_spinwave('LaMnO3.txt');
mfs = sqw_spinwave('MnFe4Si3.txt');

model = sqw_spinwave('file','powder');

model = sqw_spinwave('file.cif');

model = sqw_spinwave('cod: La Mn O3');

model = sqw_spinwave('file','edit');
model = sqw_spinwave('file','powder edit');

model = sqw_spinwave(model,'edit');

... J1=-18 ...

... J1=$J1 ...

... J1=-18 ...

... J1=($J1=-18) ...

model.J1 = -18;

p(1)=            Amplitude
p(...)=            any other parameter defined with $par symbols

• When all axes are given as vectors of same orientation and same length, the HKLE set is assumed to be list of 'events'.
• When the H K L axes are vectors of same length/orientation, the HKL locations is assumed to be a q-path in the BZ, but the E (energy) axis can be of different length/orientation to compute the Model along (q,E) path, resulting in a 2D object. Typically, the band_structure function does this, to provide the phonon dispersion curve along principal directions.
• When the H K L E axes are given as mixed orientations/length, as vectors or matrices, the Model evaluation is done on a HKLE cube, which is well suited to visualize the whole dispersion surfaces. In the example below, ql has a different orientation (it is transposed) as qh and qk, which triggers a 4D HKLE cube evaluation.

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

model(p, q, w)    % return a matrix. p can be [] to use guessed/default parameters
iData(model, p, q, w)    % return an iData

mfs = sqw_spinwave('MnFe4Si3.txt', 'edit');
d=iData(mfs, [], 0:.05:4, 0,0, 0:0.5:20);
plot(log10(squeeze(d))); % 'squeeze' removes singleton QK and QL axes

mfs = sqw_spinwave('MnFe4Si3.txt', 'powder');
q = 0.01:0.25:3; % from 0.01 to 3 Angs-1 with .25 steps
w = 0:0.5:20;
d=iData(mfs, [], q, w);
plot(d);

Phonon dispersion in perovskites

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.3-0.5 rlu and |w| < 100 meV.

S(q,w)	Description	Dimensionality	Parameters
sqw_vaks	Phonon dispersions in perovskite cubic crystals using the Vaks parameterisation	4D (HKLw)	acoustic and optical , coupling parameters, soft mode frequency. Axes in rlu.

To create the model without defined parameter values, you may use:

s = sqw_vaks;

s = sqw_vaks(p);

  p(  1)=                  At transverse acoustic slope [meV²/rlu²]
  p(  2)=                  Al longitudinal acoustic slope [meV²/rlu²]
  p(  3)=                  Aa anisotropic acoustic slope [meV²/rlu²]
  p(  4)=                  St soft mode  transverse soft optical slope [meV²/rlu²]
  p(  5)=                  Sa soft mode  anisotropic soft optical slope [meV²/rlu²]
  p(  6)=                  Vt transverse acoustic-optical coupling [meV²/rlu²]
  p(  7)=                  Va entered as: anisotropic acoustic-optical coupling [meV²/rlu²]
  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

• When all axes are given as vectors of same orientation and same length, the HKLE set is assumed to be list of 'events'.
• When the H K L axes are vectors of same length/orientation, the HKL locations is assumed to be a q-path in the BZ, but the E (energy) axis can be of different length/orientation to compute the Model along (q,E) path, resulting in a 2D object. Typically, the band_structure function does this, to provide the phonon dispersion curve along principal directions.
• When the H K L E axes are given as mixed orientations/length, as vectors or matrices, the Model evaluation is done on a HKLE cube, which is well suited to visualize the whole dispersion surfaces. In the example below, ql has a different orientation (it is transposed) as qh and qk, which triggers a 4D HKLE cube evaluation.
• When the qh,qk,ql axes span over a whole Brillouin zone, e.g. [-0.5 : 0.5], the phonon density of states (DOS) is computed, and stored in s.UserData.DOS, as well as partial DOS in s.UserData.DOS_partials (per mode).
• You may evaluate the powder dispersion using the powder(s) function.
  

	>> 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

• E. Farhi et al, EPJB, 15 (2000) pp 615-623. DOI: 10.1007/s100510051164
• A.K. Tagantsev et al, Nature Comm, 4 (2013) 2229. DOI: 10.1038/ncomms3229
• V.G. Vaks, Introduction to the Microscopic Theory of Ferroelectrics (Nauka, Moscow, 1973)

Acoustic phonon dispersion in mono-atomic cubic crystal

w(k) = E_max*sin(k.a)         along principal axes.

S(q,w)	Description	Dimensionality	Parameters
sqw_cubic_monoatomic	Phonon dispersions in a monoatomic cubic crystal using the Dynamical matrix.	4D (HKLw)	acoustic force constant ratio and scaling energy. Axes in rlu.

  p(  1)=             C_ratio C1/C2 force constant ratio first/second neighbours
  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

• When all axes are given as vectors of same orientation and same length, the HKLE set is assumed to be list of 'events'.
• When the H K L axes are vectors of same length/orientation, the HKL locations is assumed to be a q-path in the BZ, but the E (energy) axis can be of different length/orientation to compute the Model along (q,E) path, resulting in a 2D object. Typically, the band_structure function does this, to provide the phonon dispersion curve along principal directions.
• When the H K L E axes are given as mixed orientations/length, as vectors or matrices, the Model evaluation is done on a HKLE cube, which is well suited to visualize the whole dispersion surfaces. In the example below, ql has a different orientation (it is transposed) as qh and qk, which triggers a 4D HKLE cube evaluation.
• When the qh,qk,ql axes span over a whole Brillouin zone, e.g. [-0.5 : 0.5], the phonon density of states (DOS) is computed, and stored in s.UserData.DOS, as well as partial DOS in s.UserData.DOS_partials (per mode).
• You may evaluate the powder dispersion using the powder(s) function.
  

	>> s=sqw_cubic_monoatomic([ 3 3 ]); 	% create model, with C1/C2=3 E0=3
	>> 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

• E. Meisterhofer <http://lampx.tugraz.at/~hadley/ss1/phonons/scdos/sc.pdf> 2013
• Kittel, Solid State Physics, Wiley, New York, 7th ed. (1996), pp 99-103.

Phenomenological linear/quadratic dispersion

• s=sqw_linquad([H0 K0 L0 E0]);

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, half-width [meV]
  p( 16)=     DHO_Temperature   Temperature [K]
  p( 17)=          Background

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).

• When all axes are given as vectors of same orientation and same length, the HKLE set is assumed to be list of 'events'.
• When the H K L axes are vectors of same length/orientation, the HKL locations is assumed to be a q-path in the BZ, but the E (energy) axis can be of different length/orientation to compute the Model along (q,E) path, resulting in a 2D object. Typically, the band_structure function does this, to provide the phonon dispersion curve along principal directions.
• When the H K L E axes are given as mixed orientations/length, as vectors or matrices, the Model evaluation is done on a HKLE cube, which is well suited to visualize the whole dispersion surfaces. In the example below, ql has a different orientation (it is transposed) as qh and qk, which triggers a 4D HKLE cube evaluation.
• When the qh,qk,ql axes span over a whole Brillouin zone, e.g. [-0.5 : 0.5], the phonon density of states (DOS) is computed, and stored in s.UserData.DOS, as well as partial DOS in s.UserData.DOS_partials (per mode).
• You may evaluate the powder dispersion using the powder(s) function.

The Slope parameters allow to estimate the mode group velocity. Using the definition rlu = 2pi/a where a is the excitation periodicity, we can derive:

velocity_[m/s] = Slope_[meV/rlu] *a_[Angs] /2pi * 151.9

A usage example is:

• s=sqw_linquad([ 0.25 0 0 5])
• qh=linspace(0,0.5,51);qk=qh; ql=qh'; w=linspace(0.01,20,50);
• f=iData(s,[],qh,qk,ql,w);         % evaluate model into an iData
• plot3(log(f(:,1,:,:)));

Phenomenological acoustic/optic dispersion (quadratic)

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 ortho-normal reciprocal coordinate frame. The dispersion slopes apply into this frame.

The general dispersion relation has the form, with q=[ q_h, q_k, q_l ] and HKLE0 =[q_h0, q_k0, q_l0, E₀]  in [rlu and meV]::

w(q)² = E₀² + Slope1²*(q_h-q_h0)² + Slope2²*(q_k-q_k0)² + Slope3²*(q_l-q_l0)²

where the Slope parameters are specified for 3 directions. The hkl directions can be rotated by specifying the dir1 and dir2 frame vectors. The dir3 is built as dir1 x dir2.  Along given directions the dispersion is quadratic (becomes linear far from HKL0).

Such a model can be used to model an acoustic/linear dispersion (with a given slope/velocity), which has a given energy minimum E0 at a given HKL0 location, use:

• s=sqw_acoustopt([H0 K0 L0 E0]);

To model an acoustic mode, use e.g. H0=K0=L0=E0=0 and adjust the Slope parameters.
To model an optical mode, use e.g. H0=K0=L0=0 but set E0 to a finite energy (e.g. zone centre energy).
To model something else, use other parameters... ;-)

The sqw_linquad model (above) is somewhat similar, but specifies a specific HKLE location instead of the minimum of the dispersion from which the expansion is quadratic.

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

The model parameters are the following:

  p(  1)=            DC_Hdir1   Slope1 dispersion direction, 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 [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, half-width [meV]
  p( 16)=     DHO_Temperature   Temperature [K]
  p( 17)=          Background

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).

• When all axes are given as vectors of same orientation and same length, the HKLE set is assumed to be list of 'events'.
• When the H K L axes are vectors of same length/orientation, the HKL locations is assumed to be a q-path in the BZ, but the E (energy) axis can be of different length/orientation to compute the Model along (q,E) path, resulting in a 2D object. Typically, the band_structure function does this, to provide the phonon dispersion curve along principal directions.
• When the H K L E axes are given as mixed orientations/length, as vectors or matrices, the Model evaluation is done on a HKLE cube, which is well suited to visualize the whole dispersion surfaces. In the example below, ql has a different orientation (it is transposed) as qh and qk, which triggers a 4D HKLE cube evaluation.
• When the qh,qk,ql axes span over a whole Brillouin zone, e.g. [-0.5 : 0.5], the phonon density of states (DOS) is computed, and stored in s.UserData.DOS, as well as partial DOS in s.UserData.DOS_partials (per mode).
• You may evaluate the powder dispersion using the powder(s) function.

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:

velocity_[m/s] = Slope_[meV/rlu] *a_[Angs] /2pi * 151.9

A usage example is:

• s=sqw_acoustopt([ 0.25 0 0 5])
• qh=linspace(0,0.5,51);qk=qh; ql=qh'; w=linspace(0.01,20,50);
• f=iData(s,[],qh,qk,ql,w);         % evaluate model into an iData
• plot3(log(f(:,1,:,:)));

Phonon dispersion from ab-initio force estimate and dynamical matrix

The sqw_phonons model computes phonon dispersions (lattice dynamics) using the Atomic Simulation Environment (ASE) in the so-called 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, ...

S(q,w)	Description	Dimensionality	Parameters
sqw_phonons	Phonon dispersions from the Dynamical matrix, using forces estimated by ab-initio using ASE and a selection of DFT codes (EMT, GPAW, ABINIT, Elk, QuantumEspresso, VASP).	4D (HKLw)	Creation: POSCAR, CIF, PDB, ... Then, only the DHO line shape. ab-initio implies no (few) tunable parameter. Axes in rlu.

As the Phonon calculation is by itself a full story, we provide dedicated help pages:

• Models Phonons usage (technical documentation, installation)
• Models Phonons Tutorial (examples)

Model builder: defining a new model easy

>> h=ifitmakefunc;
>> h=edit(iFunc);	% equivalent

>> [p,criteria,message,output]= fits(a, h);
>> [p,criteria,message,output]= fits(h, a); % same as above

>> h=ifitmakefunc('p(1)*exp( (x-p(2))/p(3) )');
>> fits(a,  iFunc('p(1)*exp( (x-p(2))/p(3) )')); % same as above

>> h = gauss + lorz
>> h = convn(lorz, 3)
>> h = convn(gauss, lorz)
>> h = gauss .* lorz; % a 1D model
>> h = gauss  * lorz; % a 2D model

>> h=ifitmakefunc('expression');
>> 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';
>> 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, ...)
% Description
% Parameters: [fun.Parameters]
Constraint;
y=Expression;

Creating complex functions from simple functions

>> h=gauss+lorz; h.Constraint = 'p(8)=0';

Using signal convolution/correlation in new functions

>> c = fconv(a,b);		% returned convoluted object with size which is size(a)+size(b)+1
>> 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');
>> c = fconvn(a,b);	% same as above in a shorter call

• conv: convolution (calls fconv)
  
• convn: normalized convolution (calls fconvn)
  
• xcorr: cross-correlation (calls fconv in correlation mode)
  

>> a = convn(lorz, 3)		% convolution of a Lorentzian with a Gaussian of width 3
>> 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

How to write manually a model function

>> edit gauss
>> edit voigt
>> ifitmakefunc(voigt)	% pop-up a dialogue
>> edit(voigt)		% same as ifitmakefunc for iFunc objects

E. Farhi - iFit/fit models - Nov. 27, 2018 2.0.2 - back to Main iFit Page