This notebook contains a tutorial to perform MD simulations of an Octamethyl Silsesquioxane molecule with NAMD and the AMBER force field

me8t8

Table of Contents

Octa-methyl Silsesqioxane
Working with python scripts
Fitting of Methyl dihedral barrier of Octamethyl Silsesqiuoxane
Exercise: Obtaining an initial guess through a $\chi^2$ plot
View/Edit the contents on an ASCII file

(Top)

Octa-methyl Silsesqioxane (mPOSS)

The mPOSS molecule is composed of a cubic cage where silicon atoms occupy the cube vertices, and oxygen atoms are located in the cube edges (see below figure). Thus, each Si atom has a tetrahedral coordination to three O atoms and one methyl group. Methyl substitution by other chemical species makes POSS molecules highly versatile, with applications as organic solvents, polymer dispersants, catalysts, nanocomposites, diodes, and many other uses. In particular, mPOSS has found application as a coating for carbon fibers and low-dielectric films.

mPOSS molecule composed of Si (yellow), O (red), C (cyan), and H (grey) atoms. Nine different chains of consecutive O-Si-C-H covalent bonds can be constructed for each methyl group, due to the three different oxygen and hydrogen atoms that can be selected at the extremes of the chain.

mPOSS molecule non-vibrational degrees of freedom are restricted to dicrete rotational diffusion of the methyl groups (-CH3). In the AMBER force field the barrier to rotation is described by a dihedral 4-body term.

$V(\phi)=K[1+cos(3\phi)]$
Where $\phi$ is the dihedral angle defined by one of the nine combinations that can be formed with the four linked atoms O, Si, C, and H.

(Top)

Working with python scripts

All the fitting functions, settings, and tools can be carried out via a script in lieu of the graphical interface. In fact, there are complex fitting functions (like DSFinterp1DFit) that can only be run through a script.

Scripts are written in the python language. Any text editor can be used but for convenience, we will use the script editor available in MantidPlot. It has the advantage that all algorithms are automatically loaded into the namespace, thus simplifying the coding. The script editor is opened from the View menu: View $\rightarrow$ Script Window:

The script window has an upper pannel to write your code, and a lower pannel where output is printed.

Commands in the editor are executed from the Execute menu: Execute $\rightarrow$ Execute all. There's also the option to execute only the commands that are higlighted via Execute $\rightarrow$ Execute selection:

In the output panel, only two and three are printed. Notice also the green arrow indicating the last line in the script that was executed.

(Top)

Fitting of Methyl dihedral barriers of Octamethyl Silsesqiuoxane MPOSS

Following the outlined strategy, we performed molecular dynamics simulations of a single mPOSS molecule for different values of the dihedral potential barrier.

• Simulatons are 20ns long for improved resolution of $S(Q,E)$
• Simulated values for dihedral potential barrier K are 0.02, 0.03, 0.04,...,0.14 (Kcal/mol)
• $I(Q,t)$ calculated with Sassena.

The output from the Sassena calculations are files fqt_inc_run1_rms2first.h5, located in in subdirectories $USER/SHUG2015/mantid/refinement/data/K$K/T_200/, where $K is any one of the dihedral potential barriers.

The fitting proces involves the following steps:

• Loading the experimental data ($S(Q,E)$ and resolution function)
• Loading the computed $I(Q,t)$
• Transform the computed $I(Q,t)$ to a computed $S(Q,E)$ commensurable to the experimental data (energy and Q binning, convolved)
• Write up the fitting model and obtain a initial guess
• Do the fit

All these steps are coded in python script refinement.py, which can be copied to the MantidPlot script window, and executed.

A detailed explanation of this script follows:

Working directory.

Load into MantidPlot the experimental resolution file (elastic.nxs) and the experimental structure factor file exp200K.nxs. The file contains $S(Q,E)$ for a powser sample of mPOSS at T=200K. It contains 9 spectra, each at a different Q value. ($Q=$0.3, 0.5, 0.7,..,1.9) and defined in the energy domain $E=$[-0.15, 0.15]meV.

The list of values for the dihedral potential barrier $K$ for which we did simulations

We will perform a series of commands for each of the simulated $I(Q,t)$

For instance, for simulation with $K=0.02$, load Sassena output file fqt_inc_run1_rms2first.h5 into workspace incSMK0.02.

Workspace incSMK0.02_fqt.Re contains as many spectra as the experimental $S(Q,E)$, and with the same Q-values. This facilitates comparison between simulation and experiment.

What happens if our simulation has different Q-values than the experiment? This could happen if we did the simulations before the experiment. The solution is to Rebin our simulated $I(Q,t)$ along the Q-axis, so that the binning corresponds to that of the experiment.

We showcase this problem with the simulation that was done for $K=0.11$.

• Initial binning of $I(Q,t)$ is $Q=0.2, 0.3, 0.4,..., 2.0$ (the binning width is 0.1)
• Desired binning is $Q=0.3,0.5,0.7,...,1.9$ (the binning width is 0.2)

The required steps are:

• Transpose $I(Q,t)$ to $I(t,Q)$
• Execute the Rebin algorithm on $I(t,Q)$. Params=[0.2, 0.2, 2.0] indicates a bin width $w=0.2$ (the second value), and initial $Q_i=0.2+w/2=0.3$ and a final $Q_f=2.0-w/2=1.9$.
• Transpose back to $I(Q,t)$
• Rescale by half because the final binning width is twice the initial binning width

Fourier transform $I(Q,t) \rightarrow S(Q,E)$ with algorithm SassenaFFT. It will read in group workspace incSMK0.2 and append workspace incSMK0.02_sqw

The energy binning of the computed $S(Q,E)$ goes from $E_i=-2meV$ to $E_f=2meV$ with bin width $0.1\mu eV$. This binning is different than the experimental binning ($E_i=-0.15meV$, $E_f=0.15meV$, bin-width=$0.4\mu eV$), hence the rebinning step. Note that our final binnin has the same width than experiment, but over a slightly broader dynamic range. This is neccessary for the convolution step that will be performed later.

In the picture above it is shown $S(Q=1.9,E, K=0.2)$, with a prominent elastic line and a quasi-elastic broadening (notice the Log-scale).

The elastic line represents the long-time self-correlation of the hydrogen atoms. Accurately reproducing this long-time self-correlation is the weakest point in our simulations due to the disparity in simulated and experimental environment surrounding each mPOSS molecule:

• In experiment each mPOSS molecule is surrounded by neighboring molecules $\rightarrow$ slow caged diffusion of the center of mass
• In simulations the mPOSS molecule is surrounded by vacuum, and we removed global rotations and translations to focus on the internal motions of the molecule.

Our simulations cannot reproduce the diffusive motions of the molecule CoM. Thus, our elastic line will have a higher intensity than in experiments.

Remove the elastic line from the computed $S(Q,E)$. Later we will include an elastic line in our fitting model. We will treat the elastic line as an additional fitting parameter.

• The first line will return a handle to the workspace incSMK0.02_sqw
• For each spectrum, we will change the value of the elastic line (points 499 and 500) to the values of the respective neightboring points (498 and 501)

Convolve experimental resolution with computed $S(Q,E)$ incSMK0.02_sqw and store in workspace simSMK0.02. Notice that the dynamical range is now $[-0.15, 0.15]$ meV, that of the experimental resolution.

Notice the difference in order of magnitude of intensities between the simulated signal and the experimental data. We remove this disparity with a rescaling, because it has the potential to prevent the fitting from succeeding.

The initial guess, $K=0.11$ was the default value in the force-field

loop over the 9 spectra, beginning with spectrum of index=8 (corresponding to $Q=1.9$) and ending with spectrum of index=0 (corresponding to $Q=0.3$)

This is our fitting model. Instead of using the graphical interface, we can write it out as a (complex) string.

$S_{model}(Q,E) = I_{el}\cdot R(Q,E) + I_{QE}\cdot S_{sim}(Q,E) + (aE+b)$

$I_{el}\cdot R(Q,E)$ is the elastic line which we removed from the computed $S(Q,E)$

$I_{QE}\cdot S_{sim}(Q,E)$ is the interpolated structure factor via fitting function DSFinterp1DFit. The inputs for this fitting function are:

• InputWorkspaces: a list of simulated structure factors, convolved with the experimental resolution.
• ParameterValues: a list of K-values for each of the InputWorkspaces.
• LocalRegression: perform a regression for each $S(Q,E,K)$ versus $K$.
• RegressionType, RegressionWindow: regression settings

Calling the fitting algorithm Fit. Results are saved in the following workspaces:

• Fit8_Workspaces: workspace containing experimental expectrum 8, best fit, and residuals.

• Fit8_NormalisedCovarianceMatrix: correlation matrix between the fitting parameters.
• Fit8_Parameters: optimized parameters for spectrum 8.

Get a handle (ws) to the workspace optimized parameters, and retrieve optimized values for Scaling and Intensity. These optimized values will be the initial guess for the optimization of the next spectrum. This scheme is a sequential fitting of the spectra.

Print some results of the fitting to the lower panel of the fitting window. The optimal dihedral potential barrier $K$ is quite conserved along the spectra.

The value of $K$ averaged over all spectra is $\bar{K}=0.0643$. The dihedral potential barrier is actually twice this quantity.

Each methyl gropu is subject to 9 different dihedrals thus the total activation energy to methyl rotation is $E_a=9\cdot 2\cdot \bar{K}=1.16Kcal/mol$. By comparison, $E_a$ obtained from a temperature scan of $S(Q,E)$ at the BASIS beamline was $1.22Kcal/mol$.

(Top)

Obtaining an initial guess through a $\chi^2$ plot

In this exercise you will write a python script to obtain an initial guess for the fitting procedure described in Fitting of Methyl dihedral barrier of Octamethyl Silsesqiuoxane. The steps are as follows:

• Load the simulated $S(Q,E)$ and convolve with experimental resolution (the command are already writting out in python scritp of Fitting of Methyl dihedral barrier of Octamethyl Silsesqiuoxane)
• Use algorith DSFinterp to produce 20 interpolated structure factors $S(Q,E)$, each with corresponding to a different value of $K$.
• Select spectrum with index 8 of each of the 20 structure factors and fit against spectrum with index 8 of the experimental $S(Q,E)$. Extract the goodnes of fit $\chi^2$.
• Save the $\chi^2$ values to a file and plot them against $K$. The the value of $K$ at the minimum of the $\chi^2$ plot is a very good initial guess.

Load the simulated $S(Q,E)$ and convolve with experimental resolution

Open the script window in MantidPlot, and paste the relevant lines from refinement.py. Be sure to actualize value of variable rootd.

Use algorith [DSFinterp](http://archive.mantidproject.org/DSFinterp) to produce 20 interpolated structure factors

Below are the python line to call algorith DSFinterp with all needed arguments. Most of them need to be filled up.

• Workspaces needs to be passed a python list containing the names of the simulated structure after convolved with the experimental resolutio and rescaled.
• ParameterValues needs to be passed a python list containing the values of $K$ corresponding to Workspace. NOTE: the list is make up of real numbers, not strings.
• OutputWorkspaces needs to be passed a python list containing the names of the workspaces where you will store the interpolated structure factors. For instance, somethink like [interp0,...,interp19]
• TargetParameterValues are the values of $K$ corresponding to the output workspaces. For instance, the 20 values [0.03, 0.035,..,0.125].

Fit spectrum with index 8 of the simulated structure factors against experimental spectrum with index 8

We will do a fit against experiment for each of the 20 interpolated structure factors, hence we need a loop:

It's our task to write code for those commands within the loop.

We start writing a python string for the fitting model. Our fitting model is as before:

$S_{model}(Q=1.9,E) = I_{el}\cdot R(Q=1.9,E) + I_{QE}\cdot S_{sim}(Q=1.9,E) + (aE+b)$

• $R(Q=1.9,E)$ is spectrum with index 8 of workspace elastic
• $S_{sim}(Q=1.9,E)$ is spectrum with index 8 of one our interpolated structure factors.
• (aE+b) is a linear background.

We will use fitting function TabulatedFunction, which loads a particular spectrum from a workspace. Below is the string containing the model. There are two fitting functions. The first is for the elastic line and the second is for one of the interpolated structure factors. Substitute the symbols ??? with the approriate workspace names.

Next we need to call the Fit algorithm. Inspect the line of refinement.py and then write a python line calling Fit to fit each of our 20 models agains spectrum with index 8 of the experimental structure factor in the energy range [-0.13, 0.10]. Pass the following value to argument Output to store the output of the fits

Next, get a handle to the workspace containing the fit parameters. Identify the row containing the Cost function value, which is the $\chi^2$ value, and save to variable chi2. Finally, print it along with the corresponding $K$ value

These are all the commands that have to be included withing the "for i in range(20)" loop.

When you run the script, $K$ and $chi^2$ values will be printed in two columns in the lower panel of the script window:

Copy these values and save them in a file. Plot the file with your favorite Linux plotting tool. I like xmgrace which can be invoked from the terminal:

Estimate the value of $K$ for the minimum of $\chi^2(K)$

Solution: script initial_guess.py

(Top)

View/Edit the contents on an ASCII file

For those not familiar with the Linux operating system, we enumerate here a few ways to view and/or edit an ASCII file:

1 File is shown in the terminal:

2 File is shown in a separate window