Hey everybody, it's Corin, CEO and co-founder of Rowan. Today we're going to look at orbitals calculations in Rowan. And we're going to do so by looking at the molecule azulene.
Azulene is an interesting hydrocarbon because unlike most hydrocarbons, which are clear liquids or white solids, azulene is actually blue. It's found in this cool blue mushroom here—it's the compound that makes the mushroom blue. And the reason that azulene is blue is because it has these two interesting rings, a seven-membered ring here and a five-membered ring here. And if you think back to intro organic chemistry, you can remember that seven-membered rings can be aromatic if they have a positive charge in addition to three double bonds, whereas five-membered rings can be aromatic if they have a negative charge in addition to two double bonds, because then they're both 6π systems, following the 4n+2 rule.
And so what you get in azulene is actually this really strong, sort of like, delta positive character in the larger ring and delta minus character in the smaller ring. And as a side effect of this, it also happens to be blue, which is quite cool. You can imagine a good comparison here is naphthalene, which is the same number of atoms, but sort of a 6 and a 6 instead of a 7 and a 5, and naphthalene is colorless and does not have the same dipole moment. So I thought it'd be fun to look at the orbitals of azulene for this demo and sort of maybe gain some insight into this molecule.
So to submit an orbitals calculation, we go to sort of the main submission page of Rowan that you get here when you don't have a calculation selected and you click “New Orbitals Calculation.” You'll almost always want to run an orbitals calculation on an optimized structure if you want to view the orbitals for the optimized form of that structure. There are cases where you might want to view orbitals on a non-stationary or distorted structure, but in general, it's good practice to run an optimization first. So let's go over here to the azulene optimization/frequency. We can go over here to this optimization tab and watch what happened. So we're pretty close; I input this from a SMILEs string, and it looks about right. But this just gets us all the way down to the correct geometry. And then we can check the frequencies and see that, yep, there's no negative or imaginary frequencies. These are all very normal stretching modes, indicating that this is indeed a ground-state structure.
Now the reason I didn't run this on the call is that this is a DFT optimization: we have to do orbitals with DFT in Rowan, and so we need to do the optimization at DFT if we want to get the true minimum, and so this took about half an hour to run, so I didn't want to do that on the call. In contrast, the actual orbital calculation here took a lot less time. So it only took about a minute, not even a minute. And we can see here that it's resubmitted from this previous calculation. So if we lost track of where the optimization was, we could go find it this way.
Okay, so what actually happens when we run an orbitals calculation in Rowan? So we have all these tabs up here, and what this is, is it's not just orbitals, it's a wide variety of properties that let us sort of understand where the electrons are, what they're doing, what the electronic distribution and properties of this molecule is. So the idea here is a whole toolkit of closely related methods that'll let us really dive into this molecule's electronics in detail.
So we'll start here with orbitals. So this is the HOMO, the highest occupied molecular orbital. We can see here that the occupation is 2 for this, in this column, and that the energy in Hartree is about negative 0.18. If you're more used to seeing this in eV, of course we can convert to eV, and that's about what we expect for HOMO. That seems very reasonable. And now this actually illustrates what this orbital looks like. So red and blue represent different signs. We can invert this if we want to, just for visual preference. This is where the highest occupied molecular orbital is. And we see that it's sort of predominantly localized on this ring here, on the five-membered ring, which makes sense, because we remember that there's more electrons hanging out in the five-membered ring than on the seven-membered ring, and so it's not surprising that the HOMO reflects that. And we can sort of change this cutoff to make these bigger. So if we want to see this instead, this is just visual style. You know, we're drawing an isosurface, so if we use smaller values for the cutoff, the orbitals will be bigger. If we use bigger values for the cutoff, the orbitals will get much smaller.
So we can look now also—well, let's reset that because this looks a little silly—we can look at the HOMO - 1, which is different. We can look at the LUMO, which sort of starts to pick up more antibonding character and is more localized on the seven-membered ring. We can look at the LUMO + 1. And when you submit an orbitals job, you can say how many orbitals you want to see. So if we wanted to go all the way up to the LUMO + 6 or the HOMO - 6, we could do that. We just don't do it by default because often those orbitals aren't relevant, and we just really care about the frontier orbitals. Obviously, one thing you might care about here is that band gap. And so that would just be the energy of the LUMO minus the energy of the HOMO.
So the next tab over is isosurfaces. And what this lets us view is electron density and electrostatic potential. If this job had radicals—if this were unpaired electrons—we could also view the spin density isosurface. But since azulene, all the electrons are paired, it's not a radical species, we can't view the spin density because it's zero everywhere. So we can start with the electron-density isosurface. So this just tells us where the electrons are. You can see that they're in the bonds, right where you'd expect. There's some more on this side. There's slightly fewer; maybe you can imagine you see that on this side. But it pretty much just looks like bonds. And again, we can tune the cutoff to make this look bigger or smaller.
The electrostatic potential, this is sort of the positive or negative potential put out by the charged nuclei and electrons. So this tells you the force that a test charge would feel at a given point. And again, this does sort of scale, it looks a little bit like the electron density, because if you're really close to the nuclei, you'll feel a repulsive force, and as you get farther away, that attenuates. So we make the cutoff smaller, and it gets bigger. If we make the cutoff bigger, the isosurface gets smaller.
In general, the most helpful way to view this though is to actually color the electron density isosurface by the electrostatic potential. And so this is what you'll often see in figures and papers. And if we do this, what we see is it's the same shape as just the electron density isosurface, but now instead of being gray, here we have a little bit of a different cutoff, we can see the areas where there's delta positive character in blue and delta negative character in red. None of this is quite red but we can change the cutoff to sort of exaggerate it, so there that's a little more red. I mean this actually now looks exactly like what we expect. So here around the outside we have the the delta positive from these sort of acidic or you know maybe not acidic and there's a pKa but they're actually not not terrible hydrogen-bond donors, these C–H bonds. And so you can imagine these, you know, sort of sticking to a hydrogen-bond acceptor and creating some of these C–H binding interactions we often see in binding sites or crystal structures. And then on the top we have this electron cloud. You see a region of delta minus. So this is now there's a lot of electron density above and below the ring and we see that reflected in the electrostatic potential.
Something else that's a little bit more subtle here is that, you know, it doesn't come through perhaps as clearly as it would in a different system, but if we actually look, there's more red on the five-membered ring. This is sort of orange, whereas we're just yellow and green over here on the seven-membered ring. And this reflects the exact same idea that the five-membered ring is delta minus and the seven-membered ring is delta plus. So that's sort of reflected here in the electrostatic potential, which is pretty cool. We can switch to a red–blue color map if that's more visually appealing to you, where the delta positive is blue and the delta minus is red, but instead of going through the rainbow, it sort of goes through white instead, as a neutral color. It's sort of just a matter of visual style. This is more common in the literature, though.
Some other properties here: these are different atom-centered charges. So when we talk about atom-centered charge, there's no one correct way to localize a distributed electron density onto specific atoms. So there's many different schemes: there's Mulliken charge, there's Löwdin charge, there's Hirshfeld charge, CM5, RESP, CHELPG. You know, there's a lot of these different schemes. We've put two of the more common ones in here, so Mulliken and Löwdin, and in general we recommend Löwdin as Mullikin is quite basis-set sensitive. You know, all methods are basis set sensitive to some degree, but Mullikin is notoriously basis -et sensitive. But we can switch which one to visualize by clicking on the different columns. And what we can do is just look around and see what charge do different atoms have, and so we can see different numerical values. I think probably the biggest takeaway we see from this is that the hydrogens are delta plus and the carbons are delta minus, which does match exactly what we saw from the isosurface. So that's reassuring.
The bond orders tab—this is something that I actually didn't encounter for quite a while as a computational chemist. But there are ways that you can actually quantify bond order between two species. You know, we often draw different resonance structures where benzene alternates between single and double bonds. And so we sort of, you know, you'll sometimes hear that like a C–C aromatic bond has a bond order of like 1.5 on average. But in different systems that's true or untrue to varying degrees. We can actually computationally come up with guesses for what the bond order as a decimal is for a given pair of elements.
And so if we look here, for instance, at this C–C bond, the Mayer bond order is predicted to be 1.35, and Wiberg bond order is predicted to be 1.45. Both of those are reasonable. These usually match relatively closely. Here, this bond is predicted to be quite a bit longer, though, so 1.03 and 1.20. And so we might actually have a decent standing to argue that this bond here, C5 to C9, is a longer and less double C–C bond than C4 to C5. So we do indeed see this reflected in the bond distances. This one is 1.39 Å up in the top left corner here. So this is 1.39 Å versus this one is about 1.5 Å. And so you know this is just a complementary way to quantify the extent of delocalization and bonding, which can be very useful in certain circumstances.
And then our last tab here is multipole moment. This shows us the dipole moment and the quadrupole moment of this molecule. So the dipole is a vector. The traditional way to depict it is to point from the region of positive charge to the region of negative charge. You can invert that in Rowan if you want to, but that's standard. And what we see here is that we're pointing, again, from the seven-membered ring to the five-membered ring, from positive to negative, just like we expected. And we also have the magnitude, which is 0.4 Debye, and then we have the actual components of this vector.
Something that's interesting to note is that this prediction of the dipole moment is quite wrong. So the dipole moment is experimentally 1.08 Debye, whereas here we're 0.4. I'm sure we're like, you know, almost 3x too small, or almost a third of the true value for the dipole moment. It's interesting to think about why this might be. I'm not actually sure, but my guess is that this reflects sort of the self-interaction error you get with r2SCAN. So r2SCAN is, in general, a really, really good level of theory for geometries and sort of routine work. But it doesn't have the long-range Hartree–Fock exchange that you might expect from a hybrid functional. And so it's less accurate at band gaps and things like this, and I think that's reflected in the erroneous prediction of the dipole moment, so this is actually quite an interesting test case.
And if you were doing a project focusing on azulene, one of the things you might want to do right off the bat is to check that the DFT method you're using actually can reproduce the electronic distribution as quantified via the dipole moment correctly. And if it can't, then you might consider trying a different functional or reaching out to a computational expert to figure out what functional would be appropriate. (You're of course always welcome to reach out to our team at Rowan for a consult.)
And then last but not least we have the quadrupole, which is less easily visualized, unfortunately. But we can see the trace, sort of the matrix elements, we can diagonalize it. And so if you're working on a project that depends on sort of higher order multipole moments, then this might be very useful to you. Anyhow, this is how to do orbitals calculations in Rowan. I think what I've hopefully conveyed is that there's a lot of ways to get insight at a very low computational cost. So happy computing.