Hello, hello!
If anyone could explain to me how there is less than a week
left of the long month of August, I would be most grateful.
I find myself back in my SoCal home... somewhat disoriented,
excited, distraught, happy, lonely, scared, nostalgic... yearning at turns for
comfiness&coziness and for excitement&adventure.
Mostly, I feel very unprepared about the unreal reality that,
granted the benevolence of the Fates, I shall be leaving for Budapest in five days(!!).
The last two weeks or so since I last wrote were a little
crazy, and quite wonderful.
Research ended nicely- and I shall attempt to provide an
explanation of the exciting things I was doing! :]
Brace yourself, dear Reader. We are in for a trip!
So... Our simulations for LTA were failing
quite miserably, and we believe the cause (or one major cause, anyway) to be
the presence of a dipole in our unit cell structure (we'll talk more about dipoles later).
Backing up a bit. My work was with zeolites.
Zeolites are pretty crystals, which means they have repetitive structural
units-- aka unit cells. What this means is that you can tile a bunch of unit
cells together to get a zeolite-- this not only contributes to the beauty of
staring at atom-level pictures of zeolites (which I am fond of doing), but it
also means that we can save ourselves time by exploiting the symmetry of the
system. Look! A picture!
 |
| A unit cell is here boxed in red. |
In order to run our simulations, we tile together something
like 12 to 36 unit cells (pretending we actually have a 2D system, we can see 4
unit cells in the picture) because 1) we can't afford to try to model a
super large system and 2) we can use tricks to make the system seem bigger than
it actually is (to get more realistic results).
[For those more mathematically inclined, refer to periodic boundary conditions (we’re mapping onto a torus!!)]
We have essentially two versions of the code that controls
our simulations-- one which uses Ewald summation to deal with
interactions between charges, and one which does not. For the purposes of this
discussion, it doesn’t matter what Ewald does. I’m just using it to
differentiate between two different ways of running the simulation. The
important thing here is that the version of the code which uses Ewald (which is
the version we were interested in) needs the coordinates of the unit cell to be
such that the center is (0,0). The other important thing is that using Ewald
requires that the unit cell not have a dipole.
For those of us who do not remember what a dipole is, dipoles
emerge when you have an imbalance of charge. To use my professor's example from intro chem, think of
elephants pulling on a rope tied to a tree. If you have two elephants pulling
just as hard in exactly opposite directions, the tree will not fall over. If
one is pulling harder than the other, or if they’re pulling in a way that the forces
they are exerting on the tree don’t cancel out, the tree will fall over. So
having a dipole รจ
the tree falls over.
In this case, we don’t want the tree to fall over.
So we want a structure centered at zero. However, the crystal structure does not start off this way. Instead, it starts off such that the (0,0) point is at the lower left
corner.
 |
| Original, Uncentered LTA Crystal Structure |
Ewald thinks that this is not okay, and so there is a
function in the Ewald code that “centers” the structure to make sure that the
structure is the way that it likes it.
At least theoretically.
Problem: Somehow or another this function was not working
quite properly. The original crystal structure has no dipole (we know this
because when we use the non-Ewald version of our code, all is well, and when we
use Ewald, things go terribly wrong), but after centering it, it does. This
suggests that something funky is going on... and the crystal is not, in fact,
being centered properly.
So, my job for the
last few days of summer research was to understand what was going on with the centering process— with
the goal of modifying the crystal structure such that it is centered (and dipole-less) to begin
with, and undergoes no modification when put through the centering function in
the Ewald code.
Summary:
How things are--
Crystal structure (no
dipole, uncentered) --[Ewald centering function]--> Weird, dipole-carrying
structure
How I want things to be--
Modified, centered
structure –[Ewald centering function]--> Still nice and centered structure
Alright, so now the problem-solving part.
First, some terminology.
Because I can’t look at the
structure in 3D very easily, I was looking at projections onto a plane.
So what does this mean?
Imagine you’re looking at a cube. This is 3D. Now imagine you
start to push on the top of the cube and the inside of the cube kind of collapses
so you can flatten the cube into the ground, and you’re just left with….
A square,
yes? (Yes.) That square is the projection of the cube onto a plane.
Or, if we started off with a delicious donut and for some
reason decided to squash it down instead of eat it, we’d get two concentric
circles as the projection onto what I shall arbitrarily term the xy plane:
 |
| A torus (or mathematical donut) |
 |
| Projection of torus onto plane |
(Credit goes to http://mathforum.org/mathimages/index.php/Torus)
[Of course, there’s no reason why you must squish your donut down towards the floor. You could also think
of squashing it against a wall—what would that look like?]
Now that we have that, the pictures we saw earlier were the projections of the crystal structure onto the xy plane. And, again, here is the original, uncentered structure:
And here is the output of the centering function (again,
projection onto the xy plane), which then gets used in the simulation:
 |
| Centered structure. Red = oxygen atoms, blue = silicon atoms |
See how we’re now “centered” at zero?
However, there is a problem, and it is one which would clearly
give rise to a dipole. Remember, we want things to be nice and symmetric (such
that for any elephant pulling in one direction, there is another elephant
pulling just as hard in the opposite direction to balance things out)—but they
are not!
(The problem is with the O atoms, specifically, so here we're just showing the oxygens)
If we look back at the uncentered structure, which has both Si and O atoms (they're just both in the same color [we also have these "X" atoms at the center of each "ball," but don't worry about that]), we see that we have elephant oxygens along all four edges (as well as two sets along each axis)-- but once we center it, we lose two edges of elephant oxygens, and the tree comes crashing down. Clearly something is amiss.
This means it is time to understand the behavior of the centering function!
Initially, I started trying to think in 2D, and this led to
my being very confused about what was happening, such that I didn’t get the
point of how things worked. Bad idea.
Lesson learned: when trying to understand something, it’s not
stupid or simple-minded to start off with the simple model. Au contraire! Keeping things
basic allows you to see and understand the important behavior, which will help
you out when the system becomes more complicated. At any rate, it’s what
scientists seem to believe, and it was wisdom that served me well.
So, let’s focus on 1D.
First, the function:
xnew = xold
– ROUND(xold/L)*L (where L is the length of the unit cell)
So let’s pretend we’re
looking at a unit cell (in this case, let’s just look at a line segment) of
length 8, and let’s see what happens at each “quarter point” – so what happens
to 0, 2, 4, 6, and 8.
For 0:
xnew = 0 –
ROUND(0/8)*8
= 0 – ROUND(0)*8
= 0 – 0*8
= 0
For 4:
xnew = 4 –
ROUND(4/8)*8
= 4 – ROUND(0.5)*8
= 4 – 1*8
= -4
For 8:
xnew = 8 –
ROUND(8/8)*8 = 0
For 2:
xnew = 2 –
ROUND(2/8)*8 = 2
For 6:
xnew = 6 –
ROUND(6/8)*8 = -2
To illustrate what happens, we can look at the original line,
and see where the original points end up in the centered line:
Original
line:
Transformed line:
 |
| The color coding is intended to "trace" points. |
Now, I will point out the interesting things. The
points from 0 to right up to 4 (or more generally, from 0 to right up to L/2)
end up being unchanged. The points from 4 to 8 (or L/2 to L) end up becoming
the negative half of the line (well, plus 0), with L/2 becoming –L/2. It’s like
we have cut the original line at the halfway point and moved what used to be
the upper half to the end, so that the last point (L, or 8 in thi case) matches
up with zero.
The problem is, that now we have two points that get transformed
to zero, and we have no point mapped onto 4, or L/2.
To look at more pretty pictures of projections of LTA….
Below, I noted how the transformation works. The blue is the
original input structure, and the beige is the “centered” structure. Because
the unit cell has so many nice lines of symmetry, I considered what happened to
each of the four “balls” that make up the unit cell. If you don't get what's going on, don't worry about it, and just admire the pretty picture (at least, I think it's pretty). Or use the discussion about our 1D model to try to make sense of it! :]
 |
| Once again, I used color coding to trace things. The arrows are intended to provide perspective. |
And here we see that the “missing points” in the projection
are located at Lx/2 and Ly/2, which is to be expected from our discussion.
In other words... it's like instead of having a set of elephants along each edge, as well as along the axes, we now have relocated the elephants which should be along the edges corresponding to the midpoints (Lx/2, and Ly/2) to the x and y axes (respectively). This means we have two sets of oxygen elephants located along the axes... stacked on top of each other, so to speak. This image, quite properly, is rather ridiculous.
Now I’ll fast forward the story, because I’ve probably
already long lost my few readers.
(As always, feel free to ask me about anything I've mentioned. I'd LOVE to talk about it!)
Basically, we have the points at both 0 and L mapping on to
the new zero point. So what I did was figure out which points in the centered
structure which were located at zero (or rather, along an axis) originally came
from a point whose x, y, or z value was equal to L, and then moved that point in
the centered structure from 0 to L/2. [Well, actually, I moved it just a teeeeny
bit to the left of L/2, because if you’ll remember, a point at L/2 gets moved
to –L/2 (4 got moved to -4, in our 1D example above)].
To translate into elephant language, I looked at the silly stack of elephants on the axes and relocated the elephants on top to a more dignified position-- along the empty edges.
At the end of this process I had a centered, dipoleless structure (uh, in theory…). Moreover, because this structure is already centered, upon going through the centering function, nothing changes.
To translate into elephant language, I looked at the silly stack of elephants on the axes and relocated the elephants on top to a more dignified position-- along the empty edges.
At the end of this process I had a centered, dipoleless structure (uh, in theory…). Moreover, because this structure is already centered, upon going through the centering function, nothing changes.
Here is the result of centering my new centered input structure (also, just so you know, the projections onto the yz and xz planes look the same):
Yay! Things appear to be fixed!
BUT WAIT. The problem was the dipole, and now I needed to
actually calculate the dipole and make sure it was zero. I’ll spare you the
detail, but… THE STRUCTURE STILL HAD A
DIPOLE! Granted, it was a much smaller dipole than it originally had… but
still! Not cool, man!
Unfortunately, this realization happened at nearly the end of
my last day of research, and there was no real time to do more problem-solving.
What I can say is the following. If you will recall, all
these pretty pictures are projections onto a plane—aka, squished donuts. I
would argue that squished donuts are not as enjoyable as real donuts, and it is
much the same with these projections. You lose information in the flattening
process.
To get an idea of this, imagine that cube again. Now imagine
giving the side edges of the cube a jagged cut, so when you look at it straight
ahead, you see something like:
The problem is, if you squish it toward the ground, you still
just see:
This isn’t the best example, but the point is, we can imagine
that there was important information being lost in the squishing process, such
that while the projection was nicely
symmetric (and thus dipole-less), the 3D structure still had rowdy elephants
causing problems and knocking over trees.
Unfortunately, my story ends here. I have a bit more detailed
information about the remaining problem (ask if interested), but alas, at this
point my research time was over, and I needed to get ready to get on a plane to
visit a certain Russian-enthusiast close to my heart (if very far away in
spacetime). So I still don’t understand exactly what was happening in the 3D
structure that I wasn’t seeing in the projections, and I don’t understand at all
how that problem came about.
I did learn a few things, though, including:
1) I love staring at pictures of pretty crystals. !
2) Start small when trying to understand tricky things.
3) It’s okay to feel clueless and lost and not good enough. You
are good enough. Just take a deep
breath, and start grasping onto whatever you can. Eventually things will start
to make sense.
4) It is good to start small, but don’t be surprised (or
traumatized) when you take it up a level, and suddenly things are broken again.
You’ve made progress! And you’re that much better at fixing things.
I learned this while getting lost in lines of Fortran and
suspecting I was actually stupid and incompetent and while fearing that I would
never understand what the hell was going on with these zeolite things…
But I suspect that what I learned this summer will extend far
beyond the abstract world of atomistic simulations and impact the way I carry
myself in everyday life-- in the same surprising
way that playing around with code can extend beyond the virtual world and
uncover some truth about reality.
Huzzah for computational chemistry! :D
