Lessons from Computational Chemistry!

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)

Essentially, in the orange circles we have little oxygen elephants pulling on an imaginary tree at the origin. The problem is, the other two edges are missing little oxygen elephants-- so the tree is going down! Oh noes! 

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:

 x_new = x_old – ROUND(x_old/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:                          

  x_new = 0 – ROUND(0/8)*8

        = 0 – ROUND(0)*8

        = 0 – 0*8

        = 0

For 4:

  x_new = 4 – ROUND(4/8)*8

        = 4 – ROUND(0.5)*8

        = 4 – 1*8

        = -4

For 8:

 x_new = 8 – ROUND(8/8)*8 = 0

For 2:

 x_new = 2 – ROUND(2/8)*8 = 2

For 6:

 x_new = 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.

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