http://www.geo.arizona.edu/xtal/movies/NaSiOSi.html

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

Here is some Matlab from a diary file:

>> pwd ans = /export/home2/ramon/glg410matlab/Na >> ls ans = Na24.grd Na24noH.grd na nalog naorg >> 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:

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.

>> size(naorg) ans = 4900 4 >> naorg_sq = reshape(naorg',140,140); >> size(naorg_sq) ans = 140 140Luckily, 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 = 0.0683 >> x = [-4.75:xstep:4.75]; >> length(x) ans = 140 >> 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')

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')

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

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

Pretty cool, eh?

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 |

>> 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 1finds elements of a that are equal to b.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

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.

I would be curious what the plot of the difference between the 2 electron density maps looked like.Here is my analysis

Pages maintained by

Prof. Ramón Arrowsmith

Last modified November 18, 1999