diff --git a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScan.py b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScan.py
index 85a0766276f4d47446b22ed307b4da0a40ea9aa2..53c949d072e0da373505deb4a3c5832931ea69c7 100644
--- a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScan.py
+++ b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScan.py
@@ -195,8 +195,8 @@ class EnergyWindowScan(DataProcessorAlgorithm):
         delete_alg.setProperty("Workspace", self._scan_ws + '_inel_elf')
         delete_alg.execute()
 
-        x_values = mtd[self._scan_ws + '_el_eq2'].readX(0)
-        num_hist = mtd[self._scan_ws + '_el_eq2'].getNumberHistograms()
+        x_values = mtd[self._scan_ws + '_el_eq1'].readX(0)
+        num_hist = mtd[self._scan_ws + '_el_eq1'].getNumberHistograms()
         if len(x_values) < 2 and self._msdfit:
             logger.warning("Unable to perform MSDFit")
             self._msdfit = False
@@ -204,8 +204,8 @@ class EnergyWindowScan(DataProcessorAlgorithm):
         if self._msdfit:
             msdfit_prog = Progress(self, start=0.9, end=1.0, nreports=4)
             msdfit_prog.report('msdFit')
-            msd_alg = self.createChildAlgorithm("MSDFit", enableLogging=False)
-            msd_alg.setProperty("InputWorkspace", self._scan_ws + '_el_eq2')
+            msd_alg = self.createChildAlgorithm("MSDFit", enableLogging=True)
+            msd_alg.setProperty("InputWorkspace", self._scan_ws + '_el_eq1')
             msd_alg.setProperty("Xstart", x_values[0])
             msd_alg.setProperty("XEnd", x_values[len(x_values) - 1])
             msd_alg.setProperty("SpecMin", 0)
diff --git a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MSDFit.py b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MSDFit.py
index 0cabaada969696f1d0b8bd7f0189abc385c8b905..9ce70199c96d00bd176eed30a45fd05e7326d1cd 100644
--- a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MSDFit.py
+++ b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MSDFit.py
@@ -1,13 +1,14 @@
 # pylint: disable=no-init
 from __future__ import (absolute_import, division, print_function)
-from mantid.simpleapi import PlotPeakByLogValue, ConvertTableToMatrixWorkspace, SortXAxis, AppendSpectra
+from mantid.simpleapi import *
 from mantid.api import *
 from mantid.kernel import *
-from six.moves import range
+from six.moves import range  # pylint: disable=redefined-builti
 
 
 class MSDFit(DataProcessorAlgorithm):
     _output_fit_ws = None
+    _model = None
     _spec_range = None
     _x_range = None
     _input_ws = None
@@ -18,11 +19,14 @@ class MSDFit(DataProcessorAlgorithm):
         return 'Workflow\\MIDAS'
 
     def summary(self):
-        return 'Fits log(intensity) vs Q-squared to obtain the mean squared displacement.'
+        return 'Fits Intensity vs Q for 3 models to obtain the mean squared displacement.'
 
     def PyInit(self):
         self.declareProperty(MatrixWorkspaceProperty('InputWorkspace', '', direction=Direction.Input),
                              doc='Sample input workspace')
+        self.declareProperty(name='Model', defaultValue='Gauss',
+                             validator=StringListValidator(['Gauss', 'Peters', 'Yi']),
+                             doc='Model options : Gauss, Peters, Yi')
 
         self.declareProperty(name='XStart', defaultValue=0.0,
                              doc='Start of fitting range')
@@ -87,12 +91,38 @@ class MSDFit(DataProcessorAlgorithm):
 
     def PyExec(self):
         self._setup()
+        progress = Progress(self, 0.0, 0.05, 3)
+        self._original_ws = self._input_ws
 
-        # Fit line to each of the spectra
-        function = 'name=LinearBackground, A0=0, A1=0'
+        rename_alg = self.createChildAlgorithm("RenameWorkspace", enableLogging=False)
+        rename_alg.setProperty("InputWorkspace", self._input_ws)
+        rename_alg.setProperty("OutputWorkspace", self._input_ws + "_" + self._model)
+        rename_alg.execute()
+        self._input_ws = self._input_ws + "_" + self._model
         input_params = [self._input_ws + ',i%d' % i for i in range(self._spec_range[0],
                                                                    self._spec_range[1] + 1)]
+
+        # Fit line to each of the spectra
+        if self._model == 'Gauss':
+            logger.information('Model : Gaussian approximation')
+            function = 'name=MsdGauss, Height=1.0, MSD=0.1'
+            function += ',constraint=(Height>0.0, MSD>0.0)'
+            params_list = ['Height', 'MSD']
+        elif self._model == 'Peters':
+            logger.information('Model : Peters & Kneller')
+            function = 'name=MsdPeters, Height=1.0, MSD=1.0, Beta=1.0'
+            function += ',constraint=(Height>0.0, MSD>0.0, 100.0>Beta>0.3)'
+            params_list = ['Height', 'MSD', 'Beta']
+        elif self._model == 'Yi':
+            logger.information('Model : Yi et al')
+            function = 'name=MsdYi, Height=1.0, MSD=1.0, Sigma=0.1'
+            function += ',constraint=(Height>0.0, MSD>0.0, Sigma>0.0)'
+            params_list = ['Height', 'MSD', 'Sigma']
+        else:
+            raise ValueError('No Model defined')
+
         input_params = ';'.join(input_params)
+        progress.report('Sequential fit')
         PlotPeakByLogValue(Input=input_params,
                            OutputWorkspace=self._output_msd_ws,
                            Function=function,
@@ -111,58 +141,53 @@ class MSDFit(DataProcessorAlgorithm):
         rename_alg.setProperty("OutputWorkspace", self._output_param_ws)
         rename_alg.execute()
 
-        params_table = mtd[self._output_param_ws]
-
-        # MSD value should be positive, but the fit output is negative
-        msd = params_table.column('A1')
-        for i, value in enumerate(msd):
-            params_table.setCell('A1', i, value * -1)
+        progress.report('Create output files')
 
         # Create workspaces for each of the parameters
         parameter_ws_group = []
-
-        # A0 workspace
-        ws_name = self._output_msd_ws + '_A0'
-        parameter_ws_group.append(ws_name)
-        ConvertTableToMatrixWorkspace(self._output_param_ws,
-                                      OutputWorkspace=ws_name,
-                                      ColumnX='axis-1',
-                                      ColumnY='A0',
-                                      ColumnE='A0_Err',
-                                      EnableLogging=False)
-        xunit = mtd[ws_name].getAxis(0).setUnit('Label')
-        xunit.setLabel('Temperature', 'K')
-        SortXAxis(InputWorkspace=ws_name,
-                  OutputWorkspace=ws_name,
-                  EnableLogging=False)
-
-        # A1 workspace
-        ws_name = self._output_msd_ws + '_A1'
-        parameter_ws_group.append(ws_name)
-        ConvertTableToMatrixWorkspace(self._output_param_ws,
-                                      OutputWorkspace=ws_name,
-                                      ColumnX='axis-1',
-                                      ColumnY='A1',
-                                      ColumnE='A1_Err',
-                                      EnableLogging=False)
-        xunit = mtd[ws_name].getAxis(0).setUnit('Label')
+        for par in params_list:
+            ws_name = self._output_msd_ws + '_' + par
+            parameter_ws_group.append(ws_name)
+            convert_alg = self.createChildAlgorithm("ConvertTableToMatrixWorkspace", enableLogging=False)
+            convert_alg.setProperty("InputWorkspace", self._output_param_ws)
+            convert_alg.setProperty("OutputWorkspace", ws_name)
+            convert_alg.setProperty("ColumnX", 'axis-1')
+            convert_alg.setProperty("ColumnY", par)
+            convert_alg.setProperty("ColumnE", par + '_Err')
+            convert_alg.execute()
+            mtd.addOrReplace(ws_name, convert_alg.getProperty("OutputWorkspace").value)
+
+        append_alg = self.createChildAlgorithm("AppendSpectra", enableLogging=False)
+        append_alg.setProperty("InputWorkspace1", self._output_msd_ws + '_' + params_list[0])
+        append_alg.setProperty("InputWorkspace2", self._output_msd_ws + '_' + params_list[1])
+        append_alg.setProperty("ValidateInputs", False)
+        append_alg.setProperty("OutputWorkspace", self._output_msd_ws)
+        append_alg.execute()
+        mtd.addOrReplace(self._output_msd_ws, append_alg.getProperty("OutputWorkspace").value)
+        if len(params_list) > 2:
+            append_alg.setProperty("InputWorkspace1", self._output_msd_ws)
+            append_alg.setProperty("InputWorkspace2", self._output_msd_ws + '_' + params_list[2])
+            append_alg.setProperty("ValidateInputs", False)
+            append_alg.setProperty("OutputWorkspace", self._output_msd_ws)
+            append_alg.execute()
+            mtd.addOrReplace(self._output_msd_ws, append_alg.getProperty("OutputWorkspace").value)
+        for par in params_list:
+            delete_alg.setProperty("Workspace", self._output_msd_ws + '_' + par)
+            delete_alg.execute()
+
+        progress.report('Change axes')
+        # Sort ascending x
+        sort_alg = self.createChildAlgorithm("SortXAxis", enableLogging=False)
+        sort_alg.setProperty("InputWorkspace", self._output_msd_ws)
+        sort_alg.setProperty("OutputWorkspace", self._output_msd_ws)
+        sort_alg.execute()
+        mtd.addOrReplace(self._output_msd_ws, sort_alg.getProperty("OutputWorkspace").value)
+        # Create a new x axis for the Q and Q**2 workspaces
+        xunit = mtd[self._output_msd_ws].getAxis(0).setUnit('Label')
         xunit.setLabel('Temperature', 'K')
-        SortXAxis(InputWorkspace=ws_name,
-                  OutputWorkspace=ws_name,
-                  EnableLogging=False)
-
-        AppendSpectra(InputWorkspace1=self._output_msd_ws + '_A0',
-                      InputWorkspace2=self._output_msd_ws + '_A1',
-                      ValidateInputs=False,
-                      OutputWorkspace=self._output_msd_ws,
-                      EnableLogging=False)
-        delete_alg.setProperty("Workspace", self._output_msd_ws + '_A0')
-        delete_alg.execute()
-        delete_alg.setProperty("Workspace", self._output_msd_ws + '_A1')
-        delete_alg.execute()
         # Create a new vertical axis for the Q and Q**2 workspaces
-        y_axis = NumericAxis.create(2)
-        for idx in range(1):
+        y_axis = NumericAxis.create(len(params_list))
+        for idx in range(len(params_list)):
             y_axis.setValue(idx, idx)
         mtd[self._output_msd_ws].replaceAxis(1, y_axis)
 
@@ -182,6 +207,11 @@ class MSDFit(DataProcessorAlgorithm):
         copy_alg.setProperty("OutputWorkspace", self._output_fit_ws)
         copy_alg.execute()
 
+        rename_alg = self.createChildAlgorithm("RenameWorkspace", enableLogging=False)
+        rename_alg.setProperty("InputWorkspace", self._input_ws)
+        rename_alg.setProperty("OutputWorkspace", self._original_ws)
+        rename_alg.execute()
+
         self.setProperty('OutputWorkspace', self._output_msd_ws)
         self.setProperty('ParameterWorkspace', self._output_param_ws)
         self.setProperty('FitWorkspaces', self._output_fit_ws)
@@ -191,6 +221,7 @@ class MSDFit(DataProcessorAlgorithm):
         Gets algorithm properties.
         """
         self._input_ws = self.getPropertyValue('InputWorkspace')
+        self._model = self.getPropertyValue('Model')
         self._output_msd_ws = self.getPropertyValue('OutputWorkspace')
 
         self._output_param_ws = self.getPropertyValue('ParameterWorkspace')
diff --git a/Framework/PythonInterface/plugins/functions/MsdGauss.py b/Framework/PythonInterface/plugins/functions/MsdGauss.py
new file mode 100644
index 0000000000000000000000000000000000000000..dd94819fcfb271e62eece8dfa83bf63d828a3fd0
--- /dev/null
+++ b/Framework/PythonInterface/plugins/functions/MsdGauss.py
@@ -0,0 +1,67 @@
+# pylint: disable=no-init,invalid-name
+'''
+@author Spencer Howells, ISIS
+@date December 05, 2013
+
+Copyright &copy; 2007-8 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge National Laboratory & European Spallation Source
+
+This file is part of Mantid.
+
+Mantid is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3 of the License, or
+(at your option) any later version.
+
+Mantid 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 General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+File change history is stored at: <https://github.com/mantidproject/mantid>
+Code Documentation is available at: <http://doxygen.mantidproject.org>
+'''
+from __future__ import (absolute_import, division, print_function)
+
+import math
+import numpy as np
+
+from mantid.api import IFunction1D, FunctionFactory
+
+
+# For a Gaussian distribution the elastic intensity is propotional to exp(-msd*Q^2)
+# where the mean square displacement msd = <r^2>.
+
+class MsdGauss(IFunction1D):
+    def category(self):
+        return "QuasiElastic"
+
+    def init(self):
+        # Active fitting parameters
+        self.declareParameter("Height", 1.0, 'Height')
+        self.declareParameter("MSD", 0.05, 'Mean square displacement')
+
+    def function1D(self, xvals):
+        height = self.getParameterValue("Height")
+        msd = self.getParameterValue("MSD")
+
+        xvals = np.array(xvals)
+        intensity = height * np.exp(-msd * xvals**2)
+
+        return intensity
+
+    def functionDeriv1D(self, xvals, jacobian):
+        height = self.getParameterValue("Height")
+        msd = self.getParameterValue("MSD")
+
+        for i, x in enumerate(xvals):
+            e = math.exp(-msd * x**2)
+            jacobian.set(i, 0, e)
+            jacobian.set(i, 1, -(x**2) * e * height)
+            i += 1
+
+
+# Required to have Mantid recognise the new function
+FunctionFactory.subscribe(MsdGauss)
diff --git a/Framework/PythonInterface/plugins/functions/MsdPeters.py b/Framework/PythonInterface/plugins/functions/MsdPeters.py
new file mode 100644
index 0000000000000000000000000000000000000000..37bc8ac8f486d8a649d8e1853df2c218d4687aba
--- /dev/null
+++ b/Framework/PythonInterface/plugins/functions/MsdPeters.py
@@ -0,0 +1,82 @@
+# pylint: disable=no-init,invalid-name
+'''
+@author Spencer Howells, ISIS
+@date December 05, 2013
+
+Copyright &copy; 2007-8 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge National Laboratory & European Spallation Source
+
+This file is part of Mantid.
+
+Mantid is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3 of the License, or
+(at your option) any later version.
+
+Mantid 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 General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+File change history is stored at: <https://github.com/mantidproject/mantid>
+Code Documentation is available at: <http://doxygen.mantidproject.org>
+'''
+from __future__ import (absolute_import, division, print_function)
+
+import math
+import numpy as np
+
+from mantid.api import IFunction1D, FunctionFactory
+
+
+# The model of Peters & Kneller (J Chem Phys 139 165102 2013) takes into account motional heterogeneity.
+# The elastic intensity is propotional to 1/(1+msd*Q^2/6*beta)^beta
+# where the mean square displacement msd = <r^2> and beta is the paramter of the Gamma function Gamma(beta)
+
+class MsdPeters(IFunction1D):
+    def category(self):
+        return "QuasiElastic"
+
+    def init(self):
+        # Active fitting parameters
+        self.declareParameter("Height", 1.0, 'Height')
+        self.declareParameter("MSD", 0.05, 'Mean square displacement')
+        self.declareParameter("Beta", 1.0, 'beta')
+
+    def function1D(self, xvals):
+        height = self.getParameterValue("Height")
+        msd = self.getParameterValue("MSD")
+        beta = self.getParameterValue("Beta")
+
+        xvals = np.array(xvals)
+        i1 = 1.0 + (msd * xvals**2 / (6.0 * beta))
+        i2 = np.power(i1, beta)
+        intensity = height / i2
+        return intensity
+
+    def functionDeriv1D(self, xvals, jacobian):
+        height = self.getParameterValue("Height")
+        msd = self.getParameterValue("MSD")
+        beta = self.getParameterValue("Beta")
+
+        for i, x in enumerate(xvals):
+            x_sq = x**2
+            q = msd * x_sq
+            q6 = q / 6
+            a1 = 1.0 + q6 / beta
+            a = 1.0 / math.pow(a1, beta)
+            b1 = q6 * x / beta + 1
+            b = -(x_sq / 6) * math.pow(b1, (-beta - 1))
+            cqx = q6 + beta
+            c1 = math.pow((cqx / beta), -beta)
+            c2 = q6 - cqx * math.log(cqx / beta)
+            c = c1 * c2 / cqx
+            jacobian.set(i, 0, a)
+            jacobian.set(i, 1, b * height)
+            jacobian.set(i, 2, c * height)
+
+
+# Required to have Mantid recognise the new function
+FunctionFactory.subscribe(MsdPeters)
diff --git a/Framework/PythonInterface/plugins/functions/MsdYi.py b/Framework/PythonInterface/plugins/functions/MsdYi.py
new file mode 100644
index 0000000000000000000000000000000000000000..04231fe30eaa886683969ea195b601fb44b9565e
--- /dev/null
+++ b/Framework/PythonInterface/plugins/functions/MsdYi.py
@@ -0,0 +1,78 @@
+# pylint: disable=no-init,invalid-name
+'''
+@author Spencer Howells, ISIS
+@date December 05, 2013
+
+Copyright &copy; 2007-8 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge National Laboratory & European Spallation Source
+
+This file is part of Mantid.
+
+Mantid is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3 of the License, or
+(at your option) any later version.
+
+Mantid 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 General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+File change history is stored at: <https://github.com/mantidproject/mantid>
+Code Documentation is available at: <http://doxygen.mantidproject.org>
+'''
+from __future__ import (absolute_import, division, print_function)
+
+import math
+import numpy as np
+
+from mantid.api import IFunction1D, FunctionFactory
+
+
+# The model of Yi et al(J Phys Chem B 1316 5029 2012) takes into account motional heterogeneity.
+# The elastic intensity is propotional to exp(-1/6*Q^2*msd)*(1+Q^4*sigma/72)
+# where the mean square displacement msd = <r^2> and sigma^2 is the variance of the msd.
+
+class MsdYi(IFunction1D):
+    def category(self):
+        return "QuasiElastic"
+
+    def init(self):
+        # Active fitting parameters
+        self.declareParameter("Height", 1.0, 'Height')
+        self.declareParameter("MSD", 0.05, 'Mean square displacement')
+        self.declareParameter("Sigma", 1.0, 'Sigma')
+
+    def function1D(self, xvals):
+        height = self.getParameterValue("Height")
+        msd = self.getParameterValue("MSD")
+        sigma = self.getParameterValue("Sigma")
+
+        xvals = np.array(xvals)
+        i1 = np.exp(-1.0 / (6.0 * xvals**2 * msd))
+        i2 = 1.0 + (np.power(xvals, 4) * sigma / 72.0)
+        intensity = height * i1 * i2
+
+        return intensity
+
+    def functionDeriv1D(self, xvals, jacobian):
+        height = self.getParameterValue("Height")
+        msd = self.getParameterValue("MSD")
+        sigma = self.getParameterValue("Sigma")
+
+        for i, x in enumerate(xvals):
+            q = msd * x**2
+            f1 = math.exp(-1.0 / (6.0 * q))
+            df1 = f1 / (6.0 * x * q)
+            x4 = math.pow(x, 4)
+            f2 = 1.0 + (x4 * sigma / 72.0)
+            df2 = x4 / 72.0
+            jacobian.set(i, 0, f1 * f2)
+            jacobian.set(i, 1, height * df1 * f2)
+            jacobian.set(i, 2, height * f1 * df2)
+
+
+# Required to have Mantid recognise the new function
+FunctionFactory.subscribe(MsdYi)
diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScanTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScanTest.py
index f6304b866be1188ed18effcd8a9bd7f965327009..781bd8ee17d61664623c6bae550db149a0fbc2e4 100644
--- a/Framework/PythonInterface/test/python/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScanTest.py
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/WorkflowAlgorithms/EnergyWindowScanTest.py
@@ -26,9 +26,9 @@ class EnergyWindowScanTest(unittest.TestCase):
         EnergyWindowScan(InputFiles="IRS26176.RAW", Instrument='IRIS', Analyser='graphite',
                          Reflection='002', SpectraRange='3, 50', ElasticRange='-0.5, 0',
                          InelasticRange='0, 0.5', GroupingMethod='Individual', MSDFit=True)
-        scan_ws = mtd['Scan_el_eq2_0_Workspace']
-        self.assertEqual(round(scan_ws.readY(0)[0], 7), -2.4719849)
-        self.assertEqual(round(scan_ws.readY(0)[1], 7), -2.3400994)
+        scan_ws = mtd['Scan_el_eq1_Gauss_0_Workspace']
+        self.assertEqual(round(scan_ws.readY(0)[0], 7), 0.0844171)
+        self.assertEqual(round(scan_ws.readY(0)[1], 7), 0.0963181)
 
 
 if __name__ == '__main__':
diff --git a/Testing/SystemTests/tests/analysis/ISISIndirectInelastic.py b/Testing/SystemTests/tests/analysis/ISISIndirectInelastic.py
index 4b74f53b4979b3ae88db0d16ffad0bf213e5e34e..4920c2d6ec2e04a31b73161009c689b6ab3716ad 100644
--- a/Testing/SystemTests/tests/analysis/ISISIndirectInelastic.py
+++ b/Testing/SystemTests/tests/analysis/ISISIndirectInelastic.py
@@ -664,8 +664,8 @@ class ISISIndirectInelasticElwinAndMSDFit(with_metaclass(ABCMeta, ISISIndirectIn
             SaveNexusProcessed(Filename=filename,
                                InputWorkspace=ws)
 
-        eq2_file = elwin_results[1]
-        msdfit_result = MSDFit(InputWorkspace=eq2_file,
+        eq_file = elwin_results[0]
+        msdfit_result = MSDFit(InputWorkspace=eq_file,
                                XStart=self.startX,
                                XEnd=self.endX,
                                SpecMax=1)
@@ -702,8 +702,8 @@ class OSIRISElwinAndMSDFit(ISISIndirectInelasticElwinAndMSDFit):
         self.files = ['osi97935_graphite002_red.nxs',
                       'osi97936_graphite002_red.nxs']
         self.eRange = [-0.02, 0.02]
-        self.startX = 0.195082
-        self.endX = 3.202128
+        self.startX = 0.189077
+        self.endX = 1.8141
 
     def get_reference_files(self):
         return ['II.OSIRISElwinEQ1.nxs',
@@ -720,8 +720,8 @@ class IRISElwinAndMSDFit(ISISIndirectInelasticElwinAndMSDFit):
         self.files = ['irs53664_graphite002_red.nxs',
                       'irs53665_graphite002_red.nxs']
         self.eRange = [-0.02, 0.02]
-        self.startX = 0.313679
-        self.endX = 3.285377
+        self.startX = 0.441682
+        self.endX = 1.85378
 
     def get_reference_files(self):
         return ['II.IRISElwinEQ1.nxs',
diff --git a/Testing/SystemTests/tests/analysis/reference/II.IRISMSDFit.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/II.IRISMSDFit.nxs.md5
index 5d245a008dee9d1c63cdfd6980cd8f9e22438846..a50d8b3ff02dd3181c60206285cd39e8155f0192 100644
--- a/Testing/SystemTests/tests/analysis/reference/II.IRISMSDFit.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/II.IRISMSDFit.nxs.md5
@@ -1 +1 @@
-57c1e21b06c1e694564318f4e3a728a0
\ No newline at end of file
+209d8696391ae63787663f899dc05868
diff --git a/Testing/SystemTests/tests/analysis/reference/II.OSIRISMSDFit.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/II.OSIRISMSDFit.nxs.md5
index d7a7d92cb400e12d5a8181342b0d3703b25dcfc2..331bff4fc3a52bdae2906730ea5d3fdd997d62f0 100644
--- a/Testing/SystemTests/tests/analysis/reference/II.OSIRISMSDFit.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/II.OSIRISMSDFit.nxs.md5
@@ -1 +1 @@
-d7bfe3d9fd5f9eaecb607366425cf968
\ No newline at end of file
+9bf57fc31348f685b2139db56ac3ee87
diff --git a/docs/source/algorithms/MSDFit-v1.rst b/docs/source/algorithms/MSDFit-v1.rst
index c46e9536ffda7c0acd055e6bcf352411bfde1967..96907467e77dcab04ba8fc4dbfcaa37cddd6b06a 100644
--- a/docs/source/algorithms/MSDFit-v1.rst
+++ b/docs/source/algorithms/MSDFit-v1.rst
@@ -9,11 +9,14 @@
 Description
 -----------
 
-Fits :math:`log(intensity)` vs :math:`Q^{2}` with a straight line for each run
+Fits :math:`intensity` vs :math:`Q` with a straight line for each run
 to obtain the mean square displacement for a given range of runs.
 
-This algorithm operates on the QSquared workspace (*_q2*) generated by the
-:ref:`ElasticWindowMultiple <algm-ElasticWindowMultiple>` algorithm.
+This algorithm operates on the :math:`Q` workspace (*_eq*) generated by
+the :ref:`ElasticWindowMultiple <algm-ElasticWindowMultiple>` algorithm.
+
+The model used for obtaining the mean squared displacement can be
+selected. These models include 'Gaussian', 'Peters', 'Yi'.
 
 Workflow
 --------
@@ -30,22 +33,36 @@ Usage
     # Create some data that is similar to the output of ElasticWindowMultiple
     sample = CreateSampleWorkspace(Function='User Defined',
                     UserDefinedFunction='name=ExpDecay,Height=1,Lifetime=6',
-                    NumBanks=1, BankPixelWidth=1, XUnit='QSquared', XMin=0.0,
+                    NumBanks=1, BankPixelWidth=1, XUnit='momentum', XMin=0.0,
                     XMax=5.0, BinWidth=0.1)
 
-    msd, param, fit = MSDFit(InputWorkspace=sample,
-                             XStart=0.0, XEnd=5.0,
-                             SpecMin=0, SpecMax=0)
+    g_msd, g_param, g_fit = MSDFit(InputWorkspace=sample,
+                                   Model="Gauss",
+                                   XStart=0.0, XEnd=5.0,
+                                   SpecMin=0, SpecMax=0)
+
+    y_msd, y_param, y_fit = MSDFit(InputWorkspace=sample,
+                                   Model="Yi",
+                                   XStart=0.0, XEnd=5.0,
+                                   SpecMin=0, SpecMax=0)
 
-    print 'A0: ' + str(msd.readY(0))
-    print 'A1: ' + str(msd.readY(1))
+    print('Using Gauss Model')
+    print('A0: ' + str(g_msd.readY(0)))
+    print('A1: ' + str(g_msd.readY(1)))
+    print('Using Yi Model')
+    print('A0: ' + str(y_msd.readY(0)))
+    print('A1: ' + str(y_msd.readY(1)))
 
 Output:
 
 .. testoutput:: ExGeneratedDataFit
 
-    A0: [ 0.95908058]
-    A1: [ 0.11014908]
+    Using Gauss Model
+    A0: [ 0.87079958]
+    A1: [ 0.03278263]
+    Using Yi Model
+    A0: [ 0.75677983]
+    A1: [ 1.76943372]
 
 .. categories::
 
diff --git a/docs/source/diagrams/MSDFit-v1_wkflw.dot b/docs/source/diagrams/MSDFit-v1_wkflw.dot
index 7b2c9b4f30d035a8fed73a485e7502f7fa990f80..17ea6cb8c4188e41bc13bf44a44b44b0f096fdab 100644
--- a/docs/source/diagrams/MSDFit-v1_wkflw.dot
+++ b/docs/source/diagrams/MSDFit-v1_wkflw.dot
@@ -5,6 +5,7 @@ digraph MSDFit {
   subgraph params {
     $param_style
     InputWorkspace
+    Model
     XStart
     XEnd
     SpecMin
@@ -29,6 +30,7 @@ digraph MSDFit {
   }
 
   InputWorkspace              -> FitFunction
+  Model                       -> FitFunction
   SpecMin                     -> FitFunction
   SpecMax                     -> FitFunction
   FitFunction                 -> PlotPeakByLogValue
diff --git a/docs/source/fitfunctions/MsdGauss.rst b/docs/source/fitfunctions/MsdGauss.rst
new file mode 100644
index 0000000000000000000000000000000000000000..99b0fb2a38247770459741b5c410846ada8eed8a
--- /dev/null
+++ b/docs/source/fitfunctions/MsdGauss.rst
@@ -0,0 +1,26 @@
+.. _func-MsdGauss:
+
+========
+MsdGauss
+========
+
+.. index:: MsdGauss
+
+Description
+-----------
+
+The Gaussian model for Mean-squared displacement fitting is defined as:
+
+.. math:: \mbox{Height}e^{-MSD x^2}
+
+where:
+-   .. math:: `MSD` is the mean squared displacement
+-   .. math:: `Height` is the intensity scaling
+
+.. attributes::
+
+.. properties::
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/fitfunctions/MsdPeters.rst b/docs/source/fitfunctions/MsdPeters.rst
new file mode 100644
index 0000000000000000000000000000000000000000..35f8c4c47593dc37b49950fc3866ed478c9d5cf3
--- /dev/null
+++ b/docs/source/fitfunctions/MsdPeters.rst
@@ -0,0 +1,29 @@
+.. _func-MsdPeters:
+
+=========
+MsdPeters
+=========
+
+.. index:: MsdPeters
+
+Description
+-----------
+
+The Peters model for Mean-squared displacement fitting is defined as:
+
+.. math::
+
+    \frac{Height}{\power{1 + \frac{MSD x^2}{6 \beta}}{\beta}}
+
+where:
+-   .. math::`MSD` is the mean-squared displacement
+-   .. math::`\beta` is the inverse scale/rate parameter
+-   .. math::`Height` is the intensity scaling
+
+.. attributes::
+
+.. properties::
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/fitfunctions/MsdYi.rst b/docs/source/fitfunctions/MsdYi.rst
new file mode 100644
index 0000000000000000000000000000000000000000..b42be04eb7029cfde4da2ef684737067f6b090b6
--- /dev/null
+++ b/docs/source/fitfunctions/MsdYi.rst
@@ -0,0 +1,29 @@
+.. _func-MsdYi:
+
+=====
+MsdYi
+=====
+
+.. index:: MsdYi
+
+Description
+-----------
+
+The Yi model for Mean-squared displacement fitting is defined as:
+
+.. math::
+
+    Height \exp(\frac{-1}{6 x^2 MSD}) (1 + \frac{x^4 * \sigma}{72})
+
+where:
+-   .. math::`MSD` is the mean squared displacement
+-   .. math::`Sigma` is the standard deviation
+-   .. math::`Height` is the intensity scaling
+
+.. attributes::
+
+.. properties::
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/release/v3.11.0/framework.rst b/docs/source/release/v3.11.0/framework.rst
index 8795870374490750dcb0e66d9423a55831976cc2..608a455b2f428becf3d9a4f15064d2f586faba54 100644
--- a/docs/source/release/v3.11.0/framework.rst
+++ b/docs/source/release/v3.11.0/framework.rst
@@ -54,7 +54,7 @@ Improved
 - :ref:`SumSpectra <algm-SumSpectra-v1>`: Fixed a bug where a wrong fallback value would be used in case of invalid values being set for min/max worspace index, and improved input validation for those properties.
 - :ref:`LoadBBY <algm-LoadBBY-v1>`: Fixed bug where the logManager did not work with sample_name, sample_aperture and source_aperture. Also added more information regarding the sample and the selected choppers.
 - :ref:`ConvertSpectrumAxis <algm-ConvertSpectrumAxis-v2>`: Added an option to disable the sorting of the resulting axis making it useful especially for scanning workspaces. Also reduced the complexity of the operation for the default (ordered axis) case from *Nˆ2* to *N*.
-
+- :ref:`MSDFit <algm-MSDFit>` now supports model selection. Currently has the option of 3 models: MsdGauss, MsdPeters and MsdYi.
 
 Deprecated
 ##########
diff --git a/docs/source/release/v3.11.0/indirect_inelastic.rst b/docs/source/release/v3.11.0/indirect_inelastic.rst
index 2624768cbd44575e8bf96f2daffb468300b0bb95..9d99e8be05ca1c1c18464b034bac33fd63b8ca16 100644
--- a/docs/source/release/v3.11.0/indirect_inelastic.rst
+++ b/docs/source/release/v3.11.0/indirect_inelastic.rst
@@ -45,24 +45,22 @@ Bugfixes
 - Correct treatment of the resolution function: convolve sample and resolution spectra with same momentum transfer.
 - Property to pass the workspace index added to :ref:`algm-ConvolutionFitSequential`.
 
-Elwin
-~~~~~
 
-Bugfixes
---------
-- Save Result now writes to file the temperature-dependent elastic intensity normalized to the lowest temperature.
+MSDFit
+~~~~~~
 
-ConvFit
-~~~~~~~
+Improvements
+------------
+- Added model selection to MSDFit, with three current models: MsdPeters, MsdYi and MsdPeters. New models now
+  work with workspaces in Q not Q^2 (e.g. _eq workspaces 'Elastic Q')
 
-Bugfixes
---------
-- Correct treatment of the resolution function: convolve sample and resolution spectra with same momentum transfer.
-- Property to pass the workspace index added to :ref:`algm-ConvolutionFitSequential`.
 
 Jump Fit
 ~~~~~~~~
 
+General
+~~~~~~~
+
 Improvements
 ------------
 - Added a flag 'ms_enabled' (default: True) for enabling/disabling multiple scattering in vesuvio user scripts.
diff --git a/qt/scientific_interfaces/Indirect/MSDFit.cpp b/qt/scientific_interfaces/Indirect/MSDFit.cpp
index b60a2a8affd1c38ecc43477aed6a24ac40b5e5b1..a9a33357f48993c566b5c81bb1d950fac5a70194 100644
--- a/qt/scientific_interfaces/Indirect/MSDFit.cpp
+++ b/qt/scientific_interfaces/Indirect/MSDFit.cpp
@@ -73,9 +73,11 @@ void MSDFit::run() {
     return;
 
   // Set the result workspace for Python script export
+  QString model = m_uiForm.modelInput->currentText();
   QString dataName = m_uiForm.dsSampleInput->getCurrentDataName();
   m_pythonExportWsName =
-      dataName.left(dataName.lastIndexOf("_")).toStdString() + "_msd";
+      dataName.left(dataName.lastIndexOf("_")).toStdString() + "_" +
+      model.toStdString() + "_msd";
 
   QString wsName = m_uiForm.dsSampleInput->getCurrentDataName();
   double xStart = m_dblManager->value(m_properties["Start"]);
@@ -85,6 +87,7 @@ void MSDFit::run() {
 
   IAlgorithm_sptr msdAlg = AlgorithmManager::Instance().create("MSDFit");
   msdAlg->initialize();
+  msdAlg->setProperty("Model", model.toStdString());
   msdAlg->setProperty("InputWorkspace", wsName.toStdString());
   msdAlg->setProperty("XStart", xStart);
   msdAlg->setProperty("XEnd", xEnd);
diff --git a/qt/scientific_interfaces/Indirect/MSDFit.ui b/qt/scientific_interfaces/Indirect/MSDFit.ui
index 5300a0720ed8fe4c44781c7981169c77292374c3..ff46f67a2ea65a6046d7c650a2c1ee1fd1c35016 100644
--- a/qt/scientific_interfaces/Indirect/MSDFit.ui
+++ b/qt/scientific_interfaces/Indirect/MSDFit.ui
@@ -14,6 +14,34 @@
    <string>Form</string>
   </property>
   <layout class="QVBoxLayout" name="verticalLayout_5">
+   <item>
+    <widget class="QGroupBox" name="modelBox">
+     <property name="title">
+      <string>Model</string>
+     </property>
+     <layout class="QGridLayout" name="gridLayout">
+      <item row="0" column="0">
+       <widget class="QComboBox" name="modelInput">
+        <item>
+         <property name="text">
+          <string>Gauss</string>
+         </property>
+        </item>
+        <item>
+         <property name="text">
+          <string>Peters</string>
+         </property>
+        </item>
+        <item>
+         <property name="text">
+          <string>Yi</string>
+         </property>
+        </item>
+       </widget>
+      </item>
+     </layout>
+    </widget>
+   </item>
    <item>
     <widget class="QGroupBox" name="gbInput">
      <property name="title">
@@ -21,30 +49,7 @@
      </property>
      <layout class="QHBoxLayout" name="horizontalLayout_18">
       <item>
-       <widget class="MantidQt::MantidWidgets::DataSelector" name="dsSampleInput" native="true">
-        <property name="sizePolicy">
-         <sizepolicy hsizetype="Minimum" vsizetype="Fixed">
-          <horstretch>0</horstretch>
-          <verstretch>0</verstretch>
-         </sizepolicy>
-        </property>
-        <property name="workspaceSuffixes" stdset="0">
-         <stringlist>
-          <string>_eq2</string>
-         </stringlist>
-        </property>
-        <property name="fileBrowserSuffixes" stdset="0">
-         <stringlist>
-          <string>_eq2.nxs</string>
-         </stringlist>
-        </property>
-        <property name="showLoad" stdset="0">
-         <bool>false</bool>
-        </property>
-        <property name="ShowGroups" stdset="0">
-         <bool>false</bool>
-        </property>
-       </widget>
+       <widget class="MantidQt::MantidWidgets::DataSelector" name="dsSampleInput" native="true"/>
       </item>
      </layout>
     </widget>