This notebook contains a tutorial to carry out calculation of the intermediate structure factor $I(Q,t)$ with the Sassena program. It is recommended that you do the molecular dynamics tutorial for this same molecule before trying this tutorial.

Sassena calculations on a mPOSS simulation

Table of Contents

Octa-methyl Silsesqioxane
Calculation of $I(Q,t)$ and quick view with hdfview
Load Sassena output into Mantid
View/Edit the contents on an ASCII file

Section

• subsection

Examples of HTML and Markdown syntax</br>

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

In the figure above, 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 discrete 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)

Calculation of $I(Q,t)$ and quick view with hdfview

The necessary files to carry out the simulation are:

• me8t8.pdb (list of atoms and molecule name)
• hydrogens.pdb (selection of atoms over which to calculate $I(Q,t)$)
• run1_rms2first.dcd (molecular dynamics trajectory with removed global rotations and translations. 4000 frames, with frames recorded every 1ps, for a total of 4ns).
• sassena.xml (input options for the Sassena program)

If you completed the tutorial on molecular dynamics for this molecule, you should have produced file run1_rms2first.dcd yourself.

Inspect file sassena.xml (can be viewed in the web browser):

• How many Q-values are we sampling?
• For a given Q-value, how many vectors are being sampled to perform the orientational average?
• Are we calculating the coherent of the incoherent cross section? How can you tell?
• If we inspect file hydrogens.pdb (if necessary, refer to View/Edit the contents on an ASCII file), how can we tell which atoms are selected for the $I(q,t)$ calculation?

Inspect file hydrogens.pdb. Find out what atoms are being selected for the calculation of $I(Q,t)$.

All necessary files are contained in directory $HOME/SHUG2015/sassena/examples/me8t8/. Let's copy them to the scratch area:

Load the sassena module. This will make available the sassena command to actually run the calculation

Now run the calculation. We ask that you limit to four cores. It should take about 1 minute

After the calculation is done, new files are generated:

• sassena.log #several info messages regarding details of the computaion
• fqt_inc.h5 #structure factor $I(Q,t)$, in binary format.

File fqt_inc.h5 can be inspected with command hdfview. In the terminal, type:

Load file fqt_inc.h5 thought the $File \rightarrow Open$ browser.

Now double-click in the fqt block data to display the nine structure factors (one per row). Highlight the first row (row in index=0, corresponding to $I(Q=0.3A^{-1},t)$) and click in the "Line Plot" icon:

In the Line Plot Options be sure to select "Row" and then click "OK":

A simple plot of $I(Q=0.3A^{-1},t)$ will be displayed. It shows a fast decay (in less than 200ps) follow by a fluctuating plateau.

The spike at $t$~4000ps is an artifact due to poor statistics. In order to calculate $I(Q,t)$ we have to compare many pairs of frames separated in time by $t$, and there are very few pairs when $t$ becomes comparable to the simulation span. These scarcity produces an artificially high correlation, hence the spike. In general, the statistical significance degrades with increasing $t$.

(Top)

Load Sassena output into Mantid

Now we are ready to load file fqt_inc.h5 in MantidPlot.

A tutorial for Basic introduction and usage of Mantid is available here. Of special relevance is section "Loading and Displaying Data".

We start by opening MantidPlot:

The Mantid session will start. Select SNS as "Default Facility" and BASIS as "Default Instrument":

• Click on the "Manage User Directories" button.
• Click "Browse To Directory" and navigate to the location of the data files (/SNSlocal/scratch/$USER/me8t8/).
• Do the same for the default save directory.

• Click "OK".
• Click "Set".

The initial Mantid session should look like the picture below, but don't worry if some of the panels are not showing up. Panels can be brought to view from the View menu.

Mantid can plot data with the intensity normalized by the bin width, and this is the default behavior for all histogram and event data. As we are not going to use this normalization in this course we need to change the default settings. To do that, go to "View"->"Preferences..." and in the new window select "2D Plots" and untick "Normalize histograms to bin width".

As a first impression, you can think of Mantid as a collection of Algorithms and Workspaces. An algorithm is a piece of code that performs a specific task. A workspace is a (multidimensional) matrix where data are stored.

We will use algorithm "LoadSassena" :

• Type LoadSassena in the Algorithm panel and click in "Execute", the Algorithm dialog window will popup (see picture below).
• Browse to the location of the file fqt_inc.h5
• We are storing the contents of file fqt_inc.h5 into an output workspace which we name as "fqt"
• Recall that our trajectory has consecutive frames separated by 1ps. In the algorithm popup, "TimeUnit" is the separation between consecutive frames and the units are precisely picoseconds. Thus, we enter 1.0.
• Click in "Run" to do the actual loading.

Workspace "fqt" (in the figure, within the orange dashed circle) is instantiated in the workspace panel.

Let's expand the contents of fqt. If you click in the "+" sign, we see that "fqt" is actually made up of six workspaces.

• fqt_qvectors - List of generating Q-vectors, each one was orientationally averaged, each one has a different modulus (from 0.3$A^{-1}$ to 1.9$A^{-1}$).
• fqt_fq0 - This is $I(Q,t=0)$. We calculated the incoherent signal, thus $I(Q,t=0)=\sum_i |b_i^{incoh}|^2$, independent of $Q$.
• fqt_fqt.Re - Real part of $I(Q,t)$.
• fqt_fqt.Im - Imaginary part of $I(Q,t)$. For a perfect orientational average, this should be exactly zero.
• fqt_fq - $\int dt I(Q,t) = S(Q,E=0)$
• fqt_fq2 - $\int dt |I(Q,t)|^2$

Clicking in the "+" of each workspace will show a brief report of its contents. For instance, fqt_fqt.Re will inform this workspace contains 9 histograms (one per Q-value, values are 0.3$A^{-1}$, 0.5$A^{-1}$, 0.7$A^{-1}$,...,1.9$A^{-1}$) and each histogram contains 7999 bins corresponding to 7999 time points.

Let's plot the histogram corresponding to Q=0.5 and Q=0.9. First left-click in workspace fqt_fqt.Re and select "Plot Spectrum"

In the popup enter the appropriate indexes for these two Q-values,separated by a comma,then press OK.

The two histograms are depicted on the same plot.

Why does the plateau diminishes with increasing Q-value? (Hint: $I(Q,t)$ is the Fourier transform of the self-correlation function $G(r,t)$. What is the physical meaning of this correlation function?)

Notice that $I(Q,t)$ extends to negative times, but we do not have negative times in the simulation, nor in the sassena output file fqt_inc.h5. The LoadSassena algorithm assumes that $I(Q,t)$ is an even function of time, and will symmetrize the data loaded from file fqt_inc.h5. This step is justified because Newton's equations are invariant under time reversal (at least for potentials that are time-independent, as was the case of our simulation).

Exercise Plot the imaginary part of $I(Q,t)$ (workspace fqt_fqt.Im).

• How does the maximum of the imaginary part of $I(Q,t)$ compare with the maximum of the real part of $I(Q,t)$ ?
• Is the imaginary part an even function of time? Why not? (Hint: $S(Q,E)$, the Fourier transform of $I(q,t)$, is a real function)

We can inspect the actual data of workspace fqt_fqt.Re by double-clicking on the workspace. It also shows the Q-values of each histogram on the left column.

At this point, it is useful to check out the Mantid help on displaying data to be aware of other ways of plotting data, or adding a curve to and existing plot.

(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

(Top)

Section

(Top)

Subsection

(Top)

Markdown Syntax Examples

local link: link</br> remote link: http://ambermd.org font face="courier new"
$$S_{model}(Q,E)=A(Q)\cdot S_{elastic}(E) + B(Q)\cdot S_{simulation}(Q,E)\otimes S_{elastic}(E) + C(Q)+D(Q)\cdot E$$

 Quoted text 

image caption

some text