notebooks/DSFinterp1DFit/dsfinterp1Dfit_apply.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:61c0fbfd656da0101a745e69b0bd7a5f1e2631a7f05f20bfd277773d1996e736"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h1>DSFinterp1DFit tutorial<h1>\n",
      "    <h2>Fitting of Methyl dihedral barriers of Octamethyl Silsesqiuoxane</h2>\n",
      "<a id='Table of Contents'></a><h3>Table of Contents</h3>\n",
      "<a href='#Goal'>Goal -- </a></br>\n",
      "<a href='#h0.0'>Import statements -- </a></br>\n",
      "<a href='#h0.1'>Working directory -- </a></br>\n",
      "<a href='#h0.2'>Load experimental resolution -- </a></br>\n",
      "<a href='#h1.0'>Load experimental scattering law -- </a></br>\n",
      "<a href='#h2.0'>Load simulated spectra -- </a></br>\n",
      "<a href='#h3.0'>Fast Fourier transform the simulated intermediate scattering law -- </a></br>\n",
      "<a href='#h4.0'>Plot of goodness of fit versus dihedral barrier -- </a></br>\n",
      "<a href='#h5.0'>DSFinterp1DFit at all measured length scales -- </a></br>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='Goal'></a><h3>Goal</h3>\n",
      "To show the dynamical transition when the solvent is fundamentally different than water. The solvent is toluene and the solute is polyStyrene<br/>\n",
      "Toluene has a Tg~150K, while pStyrene has a Tg~370K (in the melt, for extremely long chains)<br/>\n",
      "How does the dynamical transition temperature Td changes with solvation level?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h0.0'></a><h3>Import statements</h3>\n",
      "We need to add the path to the Mantid python API, and import the Mantid algorithms  we will be using"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sys\n",
      "sys.path.append('/path/file/to/mantid/python/API') # /opt/Mantid/bin in linux systems\n",
      "from mantid.simpleapi import LoadNexus, LoadSassena, Transpose, Rebin, Scale, SassenaFFT, ConvolveWorkspaces, Fit, mtd, DSFinterp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h0.1'></a><h3>Working directory</h3>\n",
      "Define the working directory. All files should be placed here"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rootd='.'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h0.2'></a><h3>Load experimental resolution</h3>\n",
      "Input the resolution file, the elastic line.\n",
      "<center><a href=\"files/elastic.png\"><img src=\"files/elastic.png\" width=\"300\" height=\"250\" alt=\"image here\"></a></center>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "LoadNexus(Filename='{0}/elastic.nxs'.format(rootd), OutputWorkspace='elastic')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h1.0'></a><h3>Load experimental scattering law</h3>\n",
      "Load the experimental quasi-elastic spectra S(Q,E) at T=200K. It contains 9 spectra, each at different Q's (0.3, 0.5, 0.7,..,1.9), and is defined in the energy domain [-0.15, 0.15]meV\n",
      "\n",
      "<center><a href=\"files/expdata.png\"><img src=\"files/expdata.png\" width=\"700\" height=\"500\" alt=\"expdata.png\"></a></center>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      " LoadNexus(Filename='{0}/exp200K.nxs'.format(rootd), OutputWorkspace='exp200K')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h2.0'></a><h3>Load simulated spectra</h3>\n",
      "We have 14 spectra at different values of the dihedral barrier K from 0.02 to 0.15"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "parametervalues='0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15'\n",
      "for K in parametervalues.split():\n",
      "    # Load spectra into Mantid workspace. Each file contains the intermediate structure factor I(Q,t) spectra for\n",
      "    # 9 different values of Q (0.3, 0.5, 0.7,..,1.9). The names of the resulting workspaces are incSM0.03, .., incSM0.15\n",
      "    LoadSassena(Filename='{0}/K{1}/T_200/fqt_inc_run1_rms2first.h5'.format(rootd,K), TimeUnit=1.0, OutputWorkspace='incSMK{0}'.format(K))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<center><table><tr>\n",
      "<td><a href=\"files/iqt_K0.06_Q1.9_linearscale.png\"><img src=\"files/iqt_K0.06_Q1.9_linearscale.png\" width=\"300\" height=\"250\" alt=\"iqt_K0.06_Q1.9_linearscale.png\"></a> <br/>\n",
      "    <i>I(Q=1.9, E, K=0.06), Linear Scale</i></td>\n",
      "<td><a href=\"files/iqt_K0.06_Q1.9.png\"><img src=\"files/iqt_K0.06_Q1.9.png\" width=\"300\" height=\"250\" alt=\"iqt_K0.06_Q1.9.png\"></a> <br/>\n",
      "    <i>I(Q=1.9, E, K=0.06), Log-Log Scale</i></td>\n",
      "</tr></table></center>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      " What happens if one of the I(Q,t) files contains a different number of Q values? In the following special case, K=0.11, the I(Q,t) file contains spectra for values of Q (0.02, 0.03, 0.04, 0.05,...,2.0). Thus, we have to rebin the spectra in the Q-coordinate."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Transpose(InputWorkspace='incSMK0.11_fqt.Re', OutputWorkspace='incSMK0.11_fqt.Re')   # from I(Q,t) to I(t,Q)\n",
      "Rebin(InputWorkspace='incSMK0.11_fqt.Re', Params=[0.2,0.2,2.0], OutputWorkspace='incSMK0.11_fqt.Re') # Rebin in Q to (0.3, 0.5,..,1.9)\n",
      "Transpose(InputWorkspace='incSMK0.11_fqt.Re', OutputWorkspace='incSMK0.11_fqt.Re')   # from I(t,Q) back to I(Q,t)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "After rebin, the intensities I(Q,t=0) in the resulting spectra are different than for spectra at other K values. We rescale the intensities to match       "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Scale(InputWorkspace='incSMK0.11_fqt.Re', factor=0.5, Operation='Multiply', OutputWorkspace='incSMK0.11_fqt.Re')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h3.0'></a><h3>Fast Fourier transform the simulated intermediate scattering law</h3>\n",
      "Fourier transform I(Q,t) to S(Q,E). Resulting workspaces are incSMK0.03_sqw .., incSMK0.15_sqw. Notice we enter the temperature in order to apply the detailed balance correction."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "SassenaFFT(InputWorkspace='incSMK{0}'.format(K), FFTonlyRealpart=1, DetailedBalance=1, Temp=200)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rebin S(Q,E) in E to the region [-0.2, 0.2] meV, of the same order as the experimental [-0.15, 0.15] but a bit broader   "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Rebin(InputWorkspace='incSMK{0}_sqw'.format(K), Params=[-0.2,0.0004,0.2], OutputWorkspace='incSMK{0}_sqw'.format(K))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Important remark: the system simulated in one POSS molecule. thus, our simulated system does not take into account the broadening\n",
      "due to motions of the molecule center of mass that take place in the crystalline phase. As a result, the elastic line in the simulation\n",
      "is over-represented.\n",
      "\n",
      "To correct this shortcoming in the simulations we will remove the simulated elastic line from the simulated S(Q,E)\n",
      "and later include a term represented the elastic line in the fitting model.\n",
      "\n",
      "S(Q,E) for for a fixed Q (see 'SQE' below) is a list of 1000 intensity values (one intensity per energy value. The elastic line is in the middle of this list, at indexes 499 and 500. We remove these intensities just by setting these values to the adjacent ones, which are quasi-elastic in nature.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ws=mtd['incSMK{0}_sqw'.format(K)] # simulated S(Q,E)\n",
      "# Remove the elastic line for the spectra at each Q.\n",
      "for iQ in range(9):\n",
      "    Q=0.03+0.02*iQ # not neccessary but just to remind the Q-value\n",
      "    SQE=ws.dataY(ih)\n",
      "    SQE[499]=SQE[498]  \n",
      "    SQE[500]=SQE[501]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<center><table><tr>\n",
      "<td><a href=\"files/incSMK0.06_sqw.png\"><img src=\"files/incSMK0.06_sqw.png\" width=\"400\" height=\"300\" alt=\"incSMK0.06_sqw.png\"></a> <br/>\n",
      "    <i> S(Q=1.9, E, K=0.06) before removal of elastic line.<br/> Notice the Log scale on the Y-axis.</i></td>\n",
      "<td><a href=\"files/incSMK0.06_sqw_elasticLineRemoved.png\"><img src=\"files/incSMK0.06_sqw_elasticLineRemoved.png\" width=\"400\" height=\"300\" alt=\"incSMK0.06_sqw_elasticLineRemoved.png\"></a> <br/>\n",
      "    <i>S(Q=1.9, E, K=0.06) after removal of elastic line</i></td>\n",
      "</tr></table></center>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h4.0'></a><h3>Plot of goodness of fit versus dihedral barrier</h3>\n",
      "Previous to the proper fitting, we will do a Chi2 plot versus dihedral barrier K. This will give us a feeling for how the goodness of the fit changes with the value of the dihedral value K. We use algorithm DSFinterp for this purpose.\n",
      "\n",
      "We will use the structure factors  simSMK0.03, .., simSMK0.15 from the simulations to construct S(Q,E,K) with\n",
      "K running from 0.03 to 0.13 every 0.001"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "parametervalues='0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14'.split()\n",
      "inputworkspaces=['simSMK{0}'.format(K) for K in parametervalues]  # S(Q,E) from the simulations\n",
      "targetvalues=[0.03+0.001*i for i in range(101)] #target K varies from 0.03 to 0.13 every 0.001\n",
      "outputworkspaces=['out{0}'.format(K) for K in targetvalues] # interpolated S(Q,E)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "From the simulated structure factors, create interpolated structure factors for each target K"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "DSFinterp(ParameterValues=[float(K) for K in parametervalues], Workspaces=inputworkspaces, LoadErrors=0,\\\n",
      "    TargetParameters=targetvalues, OutputWorkspaces=outputworkspaces,\\\n",
      "    LocalRegression=1, RegressionWindow=4, RegressionType='quadratic')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For each target K, fit the associated structure factor S(Q=1.9, E, K)to the experimenal S(Q=1.9, E), and extract the Chi2 value from the fitting. Recall that in S(Q,E) Q goes from 0.03 to 1.9. We select the structure factor with highest Q, Q=1.9, corresponding to L~3Angstroms, to do the fitting. Among the reported Q-values, Q=1.9 corresponds to a lenght scale most sensitive to methyl rotations"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "chi2values=[]\n",
      "for K in targetvalues:\n",
      "    # Our fitting model is a elastic line at Q=1.9 plus our simulated S(Q=1.9,E) plus a linear background\n",
      "    fit_string    ='name=TabulatedFunction,Workspace=elastic,WorkspaceIndex=8,Scaling=1.0;'\n",
      "    fit_string +='name=TabulatedFunction,Workspace=out{0},WorkspaceIndex=8,Scaling=1.0;'.format(K)\n",
      "    fit_string +='name=LinearBackground,A0=0.0,A1=0.0'\n",
      "    # Algorithm Fit carries out the fitting of the previous model to the experimental S(Q=1.9,E), which is index 8 of workspace exp200K.\n",
      "    # Notice the energy range is [-0.13, 0.1] meV. Outside this energy range, the signal to noise ration of the experimental S(Q=1.9, E)\n",
      "    # is not very reliable.\n",
      "    Fit(fit_string, InputWorkspace='exp200K', WorkspaceIndex=8, StartX=-0.13, EndX=0.1, CreateOutput = 1 )\n",
      "    ws=mtd['exp200K_Parameters'] #Workspace resulting from the fitting containing the optimal value of the fitting parameters and the Chi2\n",
      "    chi2=ws.row(4)['Value']\n",
      "    print 'K=', K, ' Chi2=', chi2\n",
      "    chi2values.append(chi2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "After we retrieve the chi2 values, we don't need the interpolated structure factors:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for outws in outputworkspaces:  DeleteWorkspace(outws)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Save chi2 values to file chi2vsK.dat for inspection. Plot this file with your favorite plotting tool."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "buf='#K Chi2\\n'\n",
      "for i in range(len(targetvalues)):  buf += '{0} {1}\\n'.format(targetvalues[i], chi2values[i])\n",
      "open('/tmp/chi2vsK.dat', 'w').write(buf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<center><a href=\"files/chi2vsK.png\"><img src=\"files/chi2vsK.png\" width=\"500\" height=\"250\" alt=\"chi2vsK.png\"></a></center>\n",
      "\n",
      "The plot suggests that the dihedral value yielding a simulated S(Q=1.9,E) that compares optimally with the experimental\n",
      "S(Q=1.9, E) is Koptimal ~ 0.065. We will use this value as initial guess for fitting S(Q,E) at all Q values."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(<a href='#Table of Contents'>Top</a>)<a id='h5.0'></a><h3>DSFinterp1DFit at all measured length scales</h3>\n",
      "The model used to fit against experiment is the sum of an elastic line contribution, the simulated S(Q,E), and a linear background.\n",
      "\n",
      "<center>$S_{model}(Q,E) = a\\cdot Elastic(Q,E) + b\\cdot S(Q,E) + (cE + d)$</center>\n",
      "\n",
      "Below are the initial guess of the fitting parameters. 'Scaling' is the intensity of the elastic line, 'Intensity' is the intensity of the simulated S(Q,E), and 'K' is the dihedral barrier."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "guess={'Scaling':1.0, 'Intensity':1.0, 'K':0.065}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will use our simulated structure factors simSMK0.03, .., simSMK0.15 to search for the optimal dihedral barrier"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "parametervalues='0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14'\n",
      "inputworkspaces=' '.join( ['simSMK{0}'.format(K) for K in parametervalues.split()])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Cycle over S(Q,E) for different values of Q. We begin with the highest Q, Q=1.9, and work out way to the smallest Q, Q=0.3, doing\n",
      "a fit at S(Q,E) at each Q. The model is defined by a fit string."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for iw in [8,7,6,5,4,3,2,1,0]:\n",
      "    Q=0.3+iw*0.2\n",
      "    # Our model is a elastic line plus a  structure factor interpolated from the simulated S(Q,E), plus a linear background\n",
      "    fit_string = 'name=TabulatedFunction,Workspace=elastic,WorkspaceIndex={0},Scaling={1},constraints=(0.0001<Scaling);'.format(iw,guess['Scaling'])+\\\n",
      "        'name=DSFinterp1DFit,InputWorkspaces=\"{0}\",ParameterValues=\"{1}\",'.format(inputworkspaces,parametervalues) +\\\n",
      "        'LoadErrors=0,LocalRegression=1,RegressionType=quadratic,RegressionWindow=4,' +\\\n",
      "        'WorkspaceIndex={0},Intensity={1},TargetParameter={2},'.format(iw,guess['Intensity'],guess['K']) +\\\n",
      "        'constraints=(0.0001<Intensity);' +\\\n",
      "        'name=LinearBackground,A0=0.0,A1=0.0'\n",
      "    Fit(fit_string, InputWorkspace='exp200K', WorkspaceIndex=iw, StartX=-0.13, EndX=0.1, CreateOutput = 1 )\n",
      "    ws=mtd['exp200K_Parameters'] # Workspace resulting from the fitting containing the optimal value of the fitting parameters and the Chi2\n",
      "    Koptimal=ws.row(2)['Value']\n",
      "    print 'Q=', Q, 'L={0:5.1f}'.format( 6.2832/Q), 'K={0:6.4f}'.format(Koptimal), 'Chi2={0:3.1f}'.format(ws.row(5)['Value'])\n",
      "    # Use the optimal intensities of the elastic line and simulated quasi-elastic line as guesses for the fitting of the S(Q,E) for the next Q\n",
      "    guess['Scaling']=ws.row(0)['Value']\n",
      "    guess['Intensity']=ws.row(1)['Value']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The output of the previous print statement:\n",
      "<center><table>\n",
      "<tr> <td><b><center>Q</center></b></td> <td><b><center>L</center></b></td> <td><b><center>K*</center></b></td> <td><b><center>Chi2</center></b></td></b>  </tr>\n",
      "<tr> <td>0.3</td> <td>20.9</td> <td><b>0.0635</b></td> <td>1.8</td>  </tr>\n",
      "<tr> <td>0.5</td> <td>09.0</td> <td><b>0.0635</b></td> <td>2.4</td>  </tr>\n",
      "<tr> <td>0.7</td> <td>09.0</td> <td><b>0.0645</b></td> <td>2.0</td>  </tr>\n",
      "<tr> <td>0.9</td> <td>07.0</td> <td><b>0.0640</b></td> <td>2.7</td>  </tr>\n",
      "<tr> <td>1.1</td> <td>05.7</td> <td><b>0.0647</b></td> <td>2.4</td>  </tr>\n",
      "<tr> <td>1.3</td> <td>04.8</td> <td><b>0.0650</b></td> <td>2.5</td> </tr>\n",
      "<tr> <td>1.5</td> <td>04.2</td> <td><b>0.0651</b></td> <td>2.9</td> </tr>\n",
      "<tr> <td>1.7</td> <td>03.7</td> <td><b>0.0652</b></td> <td>3.6</td> </tr>\n",
      "<tr> <td>1.9</td> <td>03.3</td> <td><b>0.0652</b></td> <td>4.1</td> </tr>\n",
      "</table></center>\n",
      "\n",
      "\n",
      "Simulations thus predict optimal K in the range [0.0635, 0.0652], a narrow range if we take into account the wide range of length scales.\n",
      "\n",
      "<center><table><tr>\n",
      "<td><a href=\"files/fitQ0.3.png\"><img src=\"files/fitQ0.3.png\" width=\"400\" height=\"300\" alt=\"fitQ0.3.png\"></a> <br/></td>\n",
      "<td><a href=\"files/fitQ1.1.png\"><img src=\"files/fitQ1.1.png\" width=\"400\" height=\"300\" alt=\"fitQ1.1.png\"></a> <br/></td>\n",
      "<td><a href=\"files/fitQ1.9.png\"><img src=\"files/fitQ1.9.png\" width=\"400\" height=\"300\" alt=\"fitQ1.9.png\"></a> <br/></td>\n",
      "</tr></table></center>"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}