GLG410 Electron Density 3D example


As an example of Matlab's data analysis capability, I emailed Dr Robert T. Downs of the University of Arizona's Geosciences Department at the recommendation of Dr. Sharp. Professor Down's research focuses on crystal structure and bonds. Look at this animated gif of the NaSiOSi system:

I had asked:
I am looking for a dataset of x,y,parameter that
the students might operate on a bit as well as visualize with the various
contouring algorithms I expect them to learn this week.  Do you have a
convenient and simple dataset I might use to broaden our horizons?
Bob replied with this:
Here is a dataset that shows a Na atom attached to the bridging O of a
O-Si-O-SI-O linkage. The DSAA is a command that I need to import the file
into surfer. The next line, 140 140, is the grid in both x and y, so there
are 140x140 data points. The x and y are not given, only the z. The next two
lines, -4.75 4.75 is the coordinates in x and y. the next line, .001158 807.
... at eh min and max values of z.
Good luck.
I needed the good luck because I was trying to understand what the matrix scenario would be for Matlab. Here is what the original data file looks like.

I had one more question:
so the following 4 columns of data are the z values?  If so, does it go
down the first column and then to the second and so on, or is it in rows?
I get the x and y description.  What are the units of z?
With this answer:
The first 140 data points are the first row, etc.
The units are electrons/angstrom**3

Matlab data analysis

Looking at the data

First, I needed to take that funny shaped initial file, cut off the header, and transfer it to darkwing.
Here is some Matlab from a diary file:
>> pwd

ans =


>> ls

ans =


>> load na
>> load naorg
>> na_sort = sort(na);

>> plot(na_sort);
>> title('Electron density sorted')
>> xlabel('Observation#')
>> ylabel('Electrons/angstrom^3')
Here was the resulting plot:
plot 1

So I concluded that most of the data were pretty small numbers and I better be careful if I tried to contour the data because I might only see the peaks.


In order to contour, I needed to reshape the original file of electron densities into a matrix that had the same shape as that with which the observations were made:
>> size(naorg)

ans =

        4900           4
>> naorg_sq = reshape(naorg',140,140);
>> size(naorg_sq)

ans =

   140   140
Luckily, that command reshape made it easy. Note how I did it on the transposed (') version of the matrix. That was because reshape works column-wise and Bob gave me the data wrapped into rows.

Then I built the x and y vectors based upon what Bob said the distances of the observation space were:
>> xstep = (4.75- -4.75)/139

xstep =


>> x = [-4.75:xstep:4.75];
>> length(x)

ans =


>> y = x;
And then I made the contour plot:
>> contour(x,y,naorg_sq);
>> title('A few equal spaced contours')
>> xlabel('distance, angstroms')
>> ylabel('distance, angstroms')
plot 2

But you can see that I was depressed and thought maybe I made an error because there were so few contours. As I thought about it, I figured that the contours were being dominated by those few very high density observations. So, I noted that you can give the contour command a vector (in this case v) of specific contours to plot:
>> v=[0:0.01:0.5];
>> contour(x,y,naorg_sq, v);
>> title('50 lower contours')
>> xlabel('distance, angstroms')
>> ylabel('distance, angstroms')
plot 3

So that looks a bit more like the picture that Bob sent me.

Surface plots

Time to get fancy. I figured I would try to have a look at the real 3D viewing action:
>> surfl(x,y,naorg_sq)
>> shading interp
>> title('Surface plot of electron density')
>> xlabel('distance, angstroms'),ylabel('distance,angstroms'),zlabel('electrons/angstrom^3')
>> colormap pink
plot 4
Pretty cool, eh?

Clipping the data

I was not so satisfied, because obviously Bob is interested in the electron density distribution in the neighborhood, and maybe figures that the peak values are biased by actually from where the measurement was taken. We need a way to cut off those peaks.
Luckily for us, we can use MatLab's relational operators.
Relational operator Description
<less than
<=less than or equal to
>greater than
>=greater than or equal to
==equal to
~=not equal to
So let's see how they work:
>> a = 1:9, b = 9-a

a =

     1     2     3     4     5     6     7     8     9

b =

     8     7     6     5     4     3     2     1     0

>> tf = a>4

tf =

     0     0     0     0     1     1     1     1     1
finds elements lf a that are greater than 4.  zeros appear in the result 
where a <= 4 and ones appear where a>4.
>> tf = (a==b) tf = 0 0 0 0 0 0 0 0 0
finds elements of a that are equal to b.

So here is how I did it for the electron density:
>> tf = naorg_sq<0.5;
>> surf(tf)
You can see that the cylinders of 0 values are where the electron densities are greater than 0.5.

>> clip = tf.*naorg_sq;

that is the crucial step. We do a dot product (element-wise multiplication) of the two matrices and make all the electron densities that were > 0.5 = 0. But we don't want our contours or surface plot to have a bunch of holes in them, so let's keep going.

>> tf = tf-1;
>> tf = tf.*-0.5;
So there we take the tf matrix and lower everything by 1, so now the places where the electron densities were too high have values of -1 and everything else is 0. Then we multiply by -0.5, so the places with electron densities too high end up being equal to 0.5 and all else is 0. Now add the new matrix to the one we clipped moment ago and plot:
>> clip_n_clean = clip+tf;
>> surfl(x,y,clip_n_clean)
>> shading interp
>> colormap pink
>> view(285,45)
>> title('Surface plot of electron density < 0.5 e^-/ang^3')
>> xlabel('distance,angstroms'),ylabel('distance,angstroms'),zlabel('electrons/angstrom^3')
In that plot, you can see I used a surface plot with lighting, changed the shading to interpolated (works best with surfl), used the pink colormap, and changed the view position.

Here all of the plots put together onto a single figure using the subplot command.

Here is the Matlab diary for the whole thing.

What about the difference between two different datasets?

A few days later, Bob asked me:
I would be curious what the plot of the difference between the 2 electron
density maps looked like.
Here is my analysis

Laplace of Electron Density

Here is my analysis.
Pages maintained by
Prof. Ramón Arrowsmith

Last modified November 18, 1999