The following content is
provided under a Creative
Commons license.
Your support will help
MIT OpenCourseware
continue to offer high-quality,
educational resources for free.
To make a donation or
view additional materials
from hundreds of MIT courses,
visit MIT OpenCourseware
at ocw.mit.edu.
TROY VAN VOORHIS: All right.
Well, good morning, everyone.
Hope you had some at least
rest over the long weekend.
Looking forward to getting
back to talking about quantum.
So last time, which
was a whole week ago--
long time ago--
we started talking
about electronic structure.
We talked about the Born
Oppenheimer approximation
and how that gives us an
electronic Schrodinger
equation, how we look for--
what we're looking for are the
eigenvalues and eigenvectors
of that Schrodinger
equation, how
we can get lots of interesting
information out of that.
And then we broke
it down so that we
saw that there were going to
be two main knobs that we were
going to have to
turn when we were
making our approximations here.
We were going out to look at
our atomic orbital basis sets
and how we compute the energy.
And we spent the
last half of class
last time talking about
choosing basis sets.
We came up with
some understanding
of these strange names that
are written at the bottom here.
And then today, we're
going to talk about some
of these ways of
computing the energy.
And then I'm going to
flip over, and we're
going to talk about
computational chemistry tools
and how we sort of assemble
that into doing a calculation.
But again, the idea
with basis sets
was that we were
sort of questing
after this mythical complete
basis set limit where the basis
set was so big that
it was infinitely
big, that making
it bigger didn't
change your answer at all.
And we never actually get there.
But if we get close enough,
we decide that's good.
All right.
So now we're going to talk today
about computing the energy.
And thankfully,
most of the methods
for computing the energy
start from something
that I believe you are familiar
with already from 5.61.
You guys talked about the
Hartree-Fock approximation?
Yeah?
Yeah.
OK.
So-- or maybe we talked about--
did we talk about singled did
we talk about determinants?
Yes.
All right.
So there we go.
So then we talked
about a determinant.
The idea of a
determinant is that you
have a set of orbitals,
they're independent,
and then you compute
the average energy.
And so, in one
form or another, I
would write the energy of
a determinant something
like this.
I would say that
there's a piece here
which is the energy of
each electron by itself.
So these are
one-electron energies.
So for every electron in my
system, if I have n electrons,
I add up the energies of
each of those electrons.
But then, because my Hamiltonian
has electron-electron repulsion
terms, I end up with some
additional contributions
that can't be
assigned to just one
electron or another electron.
And those terms are
here in this sum.
They're pairwise so for
every mu and nu summing
over independent pairs,
I have a coulomb, term j,
and an exchange, term k.
And what these things
show me is that there
is an average repulsion
here between the electrons.
The electrons are independent.
They can't avoid each other.
But they still interact.
There's still the
electron-electron repulsion
between electrons here.
So that's an energy expression
that, in one form or another,
should be somewhat familiar,
I think, from things
you guys have done before.
The idea of Hartree-Fock
is relatively simple,
which is we have
an energy, then.
This is true for any
Slater determinant,
any set of orbitals.
So I have some orbitals here.
And the idea of Hartree-Fock
is to choose the orbitals
to minimize the energy.
This independent
particle model energy,
I want to minimize that.
I'm going to choose the
orbitals that minimize that.
And the way I'm going
to do that is I'm
going to write my orbitals--
psi mu-- as a sum
over some index, i,
some coefficients c i mu, times
my basis functions, phi i.
And I'll just emphasize
that what I have here are
functions from my
AO basis, the things
that I talked about
last time, like
the actual 1s-like, 2s-like,
2p-like, 3p-like orbitals.
And then on the other
side, what I have here
are my molecular-like orbitals.
The molecular orbitals
are linear combinations
of my atomic orbital
basis functions.
So that's the general idea
of Hartree-Fock there,
is that we're just going
to choose these orbitals
to minimize the energy.
And because the orbitals are
defined by these coefficients,
I can think of this
energy as a function
of the coefficients, c.
So I can say, all
right, that means
that I have an
independent particle model
energy which depends on
these coefficients, c.
And if I want to
minimize that energy,
I then need to
look for the places
where the derivatives
are equal to zero.
So I take the derivative of
that energy, set it equal to 0.
Now as you might suspect, if
I look at this expression,
I have an orbital here, here,
here, here, here, here, here,
here, here, here.
I have a lot of orbitals
appearing in there.
And then I'm going to have
to use the chain rule,
because every one
of those orbitals
depends on every one
of the coefficients.
So I'm going to have a
chain rule expression.
You might expect
that there's going
to be a lot of algebra involved
in taking this derivative,
and there is.
This is where I get to use my
favorite abbreviation, which
is ASA.
ASA stands for
After Some Algebra.
It means I'm not going to
work through the algebra.
There's some algebra.
There's several pages
of lecture notes
where algebra is gone through.
After you do the
algebra, you end up
with the left-hand equation
reducing to an equation that
looks like this.
So there is a matrix,
which I will call h,
and then when I dot that
matrix into a vector, c
mu, the coefficient vector,
I get some energy, e,
mu, times the
coefficient vector again.
So that takes a lot of
algebra to get to that,
but you can get an equation
that looks like this.
And that's encouraging,
because this is familiar to me.
This is an eigenvalue equation.
And that's nice because I
remember eigenvalue equations
from--
I mean, the Schrodinger equation
is an eigenvalue equation.
So it's nice that when
I do this minimization,
I get back an
eigenvalue equation.
That's nice.
The one complication is
that the Hamiltonian here
depends on the coefficient, c.
So it's some
effective Hamiltonian.
I'm intentionally
collecting terms
so that it looks like
a Schrodinger equation,
but it's a little bit funny
because the Hamiltonian here
depends on the coefficients
that solve the equation.
So this is sort of an
eigenvalue equation.
And so what we're seeing
here is that H depends
on the wave function somehow.
It seems like, depending on
what the wave function is,
I have different Hamiltonians.
And to understand
why this is the case,
I want us to think
about the potential, v,
that an electron feels in
an atom or in a molecule.
And there's generally going
to be two pieces here.
And they're actually--
I can color-code
them with the terms
that exist in the independent
particle expression.
So there is the
one-electron terms,
the terms in the energy that
just come from one electron.
And those are going to be
things like, well, I might have
a nucleus of charge plus z.
And I might have my
electron out here, the one
that I'm interested in.
And it's going to have some
attraction to that nucleus.
So that just involves
one electron.
It involves a nucleus,
but just one electron.
So there's some
nuclear potential
that each electron
is going to feel.
But then there's the
term that correspond
to these average
repulsion terms.
So there are other
electrons in my system.
And each of those electrons
interacts with the electron
that I'm looking at.
So there's going to be some
electron-electron repulsion
term here.
And that electron-electron
repulsion
depends on where the
other electrons are.
And those other electrons
have some wave function.
They're kind of spread out
here, all over the place.
And depending on how spread
out those other electrons are,
I might have more or
less electron repulsion.
So if they're very, very
concentrated in the region
that I'm looking at,
the electron repulsion
is going to be very large.
If they're very spread
out, and very diffuse,
they're not going to have
very much electron repulsion.
But overall, the
electron-electron repulsion
depends on the wave
functions that I've chosen.
So the potential here is called
a mean field, or average field.
And it's because it's
the potential that comes
from the average repulsion.
It's the effective potential due
to the averaging, the smearing
out, over all the
other electrons.
So what we do have, then,
is that the Hamiltonian
does into indeed depend
on the wave functions,
because the wave
functions determine
that average potential.
But then if I know
the Hamiltonian,
I can solve for
the wave functions,
because the wave functions
are the eigenfunctions
of the Hamiltonian.
And so I have this sort of
chicken and egg situation,
that the Hamiltonian
defines the wave functions,
and the wave functions
define the Hamiltonian.
And the solution
to this is to do
what are called self-consistent
field iterations.
And you can more or
less think about these
as the simple process of,
guess some wave functions,
build the Hamiltonian, get
some new wave functions,
build a new Hamiltonian,
get some new wave functions,
build a new
Hamiltonian, and keep
going until everything
becomes consistent.
And the idea, then,
is that the wave
function we're looking
for in Hartree-Fock
are the eigenfunctions of H of
psi Hf, which we just defined
to be equal to H Hartree-Fock.
So there's some
Hartree-Fock Hamiltonian.
And the Hartree-Fock
wave functions
are the eigenfunctions
of that Hamiltonian.
Does that makes sense?
Any questions about that?
All right.
So Hartree-Fock is
the starting point.
So we have that as
an approximation.
It's a useful starting point.
In many cases, Hartree-Fock
is not good enough.
And so the real thing that
we're trying to figure out
is, how do we make
something that
uses what we've learned from
Hartree-Fock about orbitals,
and energies, and
average repulsion-- how
do we make something that's
like Hartree-Fock, but better?
So one way to do
this is something
that you've spent a lot of
time on this semester, which
is perturbation theory.
You have something that's
pretty good as a starting point,
and you want to make it
better, perturbation theory
is the natural thing
to turn to do that.
And so the way we apply
perturbation theory
in Hartree-Fock is
we say, well, we've
got this Hartree-Fock
Hamiltonian, which
we can solve, and then we've
got the actual Hamiltonian, H,
which we can't solve.
And so what we can do is
then we say, aha, well
that means that if I
just take the difference
between the
Hamiltonian I can solve
and the Hamiltonian
I can't solve,
and I treat that
as a perturbation,
then I can do a well-defined
perturbation theory where
I know the solutions
to H0, and then I
treat that additional
term as a perturbation.
And so then if I
do that, I end up
with a perturbation
expansion where
I can go to first, second,
third, fourth, fifth order,
and that additional term that's
the difference between the two.
And if I look at the first
two terms added together,
they give me the
Hartree-Fock energy.
So if I'm the one making up
a perturbation expansion,
I will always choose it
so that the first order
correction is 0.
So that if otherwise
I have chosen,
I will just re-engineer my
perturbation so that's true.
So up through first order, I
just get what I started with.
I get Hartree-Fock.
But then I have these
additional terms.
And so I can go to second order.
And I will get a
method that I call--
it's called MP2.
So Hartree-Fock is named
after Hartree and Fock.
MP stands for
Moller and Plesset,
the names of the two people
who wrote the paper in 1938,
about this particular thing.
And the 2 is second order.
And then I could go on, and I
could do MP3, MP4, dot dot dot
dot dot.
So in general, this is
MP, and some number, n,
telling me the order I go
to in perturbation theory.
So you could do this.
In principle, this gets
better and better answers.
However, in practice,
about 20 years ago,
people discovered something
that was rather disturbing.
So if I take a look at the
energy minus the [INAUDIBLE]..
So let's say I have a case
where I know the exact energy.
And I look at the energy of
these different approximate
methods for, say, one particular
simple thing, like a neon atom.
So I take a neon
atom, and I look
at the energies of these
MPn approximations,
and I look at them as a function
of n, so 1, 2, 3, 4, 5, 6, 7,
8.
And I'll say that for neon--
neon is something for which
Hartree-Fock is already
pretty good.
So the difference
between Hartree-Fock
and the exact answer
is not too bad.
So Hartree-Fock
might be, let's say,
there on this particular scale.
And then if you do
MP2, MP2 is better.
So my MP2 is maybe there.
MP3 is also a little
bit better than that,
but it overshoots a little bit.
MP4 corrects back in
the opposite direction.
And then we've got 5, 6, 7, 8.
And you can-- if I
connect the dots,
you can see what's
starting to happen here.
So that's a series.
So mathematically,
I've got something
that there is a value
for everything in n.
That's a series.
And the name for
a series like this
is one that is not converging.
So there's an answer which
is the infinite order answer,
but these partial sums
are not converging
to the infinite answer.
And very often this turns out
to be the case with perturbation
theory.
There can be
perturbation theories
that are very useful at
low orders for describing
qualitatively how things shift.
But if you go out
to infinite order,
they are series that, in
principle, you can prove,
give the correct, exact
answer, because we've
derived it as such.
But just because a series in
principle sums to something
doesn't mean that any
finite number of terms
will give a convergent answer,
or give smaller errors.
And people find that
for perturbation theory
for electronic structure
theory, this is--
not all of the
time, but very often
the case, that the perturbation
expansion doesn't converge.
And so this is one
of those places where
you make a nice empirical
compromise, and you say,
well, gee, every term--
it doesn't necessarily
get better after MP2.
MP2 usually improves
on Hartree-Fock,
but higher order might
be better or worse.
And certainly, higher
order is more expensive.
So MP2 is a good place
to get off the elevator,
because it's just going to
get more expensive, and not
necessarily more accurate.
So this is a good place to stop.
And so you'll see a
lot of calculations
out there that are
variations on MP2,
because it's the
least expensive,
and it's somewhat more accurate.
OK so now we're going
to take a pause,
and we're going to talk--
and this is where I need
audience participation.
So we're going to discuss
two different things
that you might want to know.
We've now got two methods,
Hartree-Fock and MP2.
And there are two things
that you might want to know.
One is, how accurate will
one of these calculations be?
If I do it, will it give
me the right answer?
And the other one is,
how long will it take?
In other words, will
it finish by Friday
when the P set is due?
Will it finish by the
end of the semester?
Will it finish by the
end of my time at MIT?
What time scale
are we looking at?
So we'll start off
talking about properties.
So people have done many,
many, many calculations
using many different
electronic structure methods,
so that now there's sort of
an empirical rule of thumb
for just about any property
you want to know about.
There's a rule of thumb for how
accurate one of these methods
should be.
So what's a property that
you might want to know?
This is where the audience
participation comes in.
I will only give you
information about properties
that you say are interesting.
AUDIENCE: Time and
space complexity.
TROY VAN VOORHIS: Time
and space complexity.
All right.
Could you be more specific.
What type of time
and space complexity?
AUDIENCE: How well--
I mean, for a size of
molecule, how long it takes--
TROY VAN VOORHIS: Ah, OK.
AUDIENCE: --and also
how much space it takes.
TROY VAN VOORHIS: Yes, so
that will be the second thing
that we will talk about.
Those are the two
variables we'll
talk about, in terms of
how long and how fast.
I'm interested here,
where I'm saying, like,
molecular properties--
things you
might want to know about
a molecule, or a reaction,
or something like that.
AUDIENCE: Like the energy of a
molecule in transition state.
TROY VAN VOORHIS: OK, so that
would be sort of activation
energy, right/ OK,
what's something else?
We only want to know
about activation?
Yeah.
AUDIENCE: Charge distribution.
TROY VAN VOORHIS: Right.
OK, so I'll give one
variable that probes that,
which is a dipole
moment, for example.
You might want to
know a dipole moment.
AUDIENCE: What the MOs
look like [INAUDIBLE]
TROY VAN VOORHIS: OK, yes.
So that is something
we want to know about,
but it's a qualitative one.
So what I'm going to give is I'm
going to give you things where
I can say, oh, for
a dipole moment,
it's within this many
debye, or this much percent.
So MOs are qualitative, and we
can get that, but I can't say,
well, this MO is 75% right.
AUDIENCE: Average distance
of an electron from an atom.
TROY VAN VOORHIS: OK, so
that would be like size.
AUDIENCE: Yeah.
TROY VAN VOORHIS: Yeah, yeah.
AUDIENCE: [INAUDIBLE].
TROY VAN VOORHIS: What was that?
AUDIENCE: Color.
TROY VAN VOORHIS: Yeah.
I mean, so yeah, that's
about an absorption--
like an electronic
absorption spectrum.
Other things people
might want to know?
AUDIENCE: Bond length.
TROY VAN VOORHIS: Bond length.
There you go.
Anything else?
We've got a good list.
OK, we'll use that as our list.
OK, so I will start
off with bond lengths,
because that's the
one where Hartree-Fock
does the best of these things.
So it makes you feel
like it's encouraging.
So bond lengths--
Hartree-Fock usually
predicts bond lengths to
be a little bit too short.
So the molecules are a
little bit too compact.
But it's not too bad.
They're usually too short by
about 1%, which is not too bad.
It gets 99% percent of
the bond length right.
MP2 does better, so
that it doesn't actually
have a systematic error.
It typically gets
bond lengths correct,
plus or minus one picometer.
So a picometer is 10 to
the minus 12th meters.
So 0.01 angstroms,
since a bond length
is usually about an angstrom.
So we'll go from that to size.
So for the size of an
atom or a molecule,
this has to do with the Van
der Waal's radius, basically.
And actually, Hartree-Fock
does a fairly good job
of describing Van
der Waal's radii.
I would say that it more or
less goes with the bond length
prescription, so that the size
is a little bit too small.
So it might be minus
2% to 3% to small.
And then there's very few
people that actually do sizes
with MP2, for various reasons.
But I would say that it's--
I'll put "accurate" there.
Because basically the size of
an atom or a molecule is fuzzy.
So I can talk about the radius,
and it's a little bit fuzzy.
The errors in MP2
would be smaller
than whatever fuzziness I could
have in my definition of size.
For activation energies,
Hartree-Fock is pretty poor,
so the activation barriers here
are too high, by 30% to 50%.
So the barrier heights are
very high in Hartree-Fock.
MP2 has barrier
heights that are still
too high, but only by
about 10% on average.
Dipole moments in
Hartree-Fock-- there
are no systematic
studies on that.
I will just say
that they are bad.
They are bad enough that
nobody does a systematic study.
But for MP2, they're quite
a bit better, and typically
accurate say, 0.1 debye.
And then for excitation
absorption energy--
so this is going to be the
difference between the lowest
electronic state and the next
highest electronic state--
so you can get these things
through various means.
The typical absorption-- and
also, I'll say that doing that
involves an extension beyond
what I've talked about
so far, because I've been
focused on the ground state.
But if you use the logical
extensions of these things,
the absorption energies
in Hartree-Fock
tend to be too big, and too
big by, say, half an eV or so.
And then for MP2, they're more
accurate, but not that much.
So plus or minus, say, 0.3 eV.
Electronic inside states tend
to be fairly difficult to get.
But we see that MP2 is doing
better on all accounts,
than Hartree-Fock.
So it's improving.
That's good.
But then, all right, what
is this going to cost me?
How long is this going to take?
Well, there's two different
ways that we can measure
computational complexity.
One is by the amount of
storage that's required.
It's the amount of information
you have to store in order
to do the calculation.
The other way is by how
long it's going to take,
how many operations the computer
has to do to solve the problem.
And so for Hartree-Fock, the
main memory piece that we need
is actually the Hamiltonian.
So we need that
Hamiltonian matrix.
We've got to store that in order
to compute its eigenvalues.
And it is an end-- it's a number
of basis functions by number
of basis functions.
And so then now we have
to say, all right--
well, first of all, we have
to figure out how much RAM.
So when I'm figuring out
how much RAM is required--
does anybody know how much
RAM they have on their laptop
or desktop computer?
Nobody knows.
16 gigs.
All right, that's
what's on my mine, too.
16 gigabytes, which is 16 times
10 to the 9th bytes, roughly.
So now we have to
say, all right,
for storing n basis
times n basis numbers,
how much does that take?
Well, let's say that
we have a atoms.
So we've got-- the number
of atoms we have is a.
Then the number of basis
functions we have-- last time,
we figured out that
a DZP basis, which
is a decent-size basis
for carbon had around 15.
It was 14 basis
functions per carbon.
I'm just going round to 15
because it's easier math.
So that means that the number
of basis functions is 15,
roughly, times the
number of atoms.
So there's 15a times
15a things that I need
to store to get this matrix.
And that means I need about
225 a squared numbers.
So these are how many
numbers I have to store.
And then I have to
know that each number--
I'm usually storing these
in double precision.
So each number requires eight
bytes in double precision.
So that means that I have
something like 1,700 a squared
bytes to store that object.
So now that says, all
right, if I've only
got 16 times 10 to the
9 bytes of storage space
on my computer, and I need 1,700
a squared bytes for a molecule
with a atoms, that just
gives me a natural limit
on the number of atoms.
If you back it out, that implies
that a is less than or about
equal to 3,000, which
is a big number.
So your molecule could have
up to about 3,000 atoms in it,
and this calculation would run.
And that rough calculation
turns out to be about right.
Now let's take a look at
the CPU time required.
So the CPU time, we're going
to measure this in hours,
because that's my
human unit of time.
So one hour is 3,600 seconds.
And then what I need
to know is well,
what actually takes
the computer time?
What takes the computer time is
doing mathematical operations--
so add, subtract, multiply,
divide, if, things like that.
And usually, we measure
the speed of a CPU
in terms of the
number of operations
you can do per second.
And does anybody know
the order of magnitude
of operations per second the
CPU on your computer can do?
It's about a billion--
some number of billion
operations per second,
depending on how recent
a model you have,
and whether you play
video games or not,
you will do more
or less than that.
So that means that our computer,
in one hour, can do about one--
and I'm going to round
up, and say you've
got a really good one.
So we're going to say
that your thing can do 10
billion operations per second.
If you've got multiple cores,
they can each go independently.
So 1 times 10 to the 10
operations per second.
And that means that you
can do about 3 times 10
to 13 operations in an hour.
I'm just going to set the
hour as my patience limit.
I don't want to wait
longer than an hour.
And so, then, for Hartree-Fock,
the rate-determining step
is computing the eigenvalues
of the Hamiltonian.
And that requires
n cubed operations,
where n is the
dimension of the matrix.
So if our matrix n basis on
a side, it requires n cubed--
n basis functions
cubed operations.
And so then, again, using
our translation of that
into atoms, that's 15 times
the number of atoms cubed.
Which, if I did my math
right back in my office,
that's about 3,000 a cubed.
And then if I back that out--
again, if I compare that to
the number of operations,
that gives me a
limit on the number
of atoms that I can
handle in an hour,
and that turns out
to be about 1,000.
So similar sizes.
So I run out of storage
about the same time
as I run out of patience.
If I was willing to be a
little bit more patient,
I might be able
to do a few more.
But order of magnitude,
we're at 1,000.
So I'll summarize.
What we found here is that
the storage requirements
were something like a
squared of order a squared.
CPU time was of order a cubed.
The maximum feasible
number of basis functions
that I could do
here was something
like-- so if I have 1,000
atoms, that was my limit.
Something like 15,000 would
be the maximum feasible
number of basis functions.
And maximum feasible number
of atoms was about 1,000.
So we won't go through
the same exercise for MP2
in such detail.
I will just tell you that
everything is worse for MP2.
Everything takes longer,
it takes more storage,
everything's worse.
So it requires order
a cubed storage.
It requires order a
to the fifth CPU time.
The maximum feasible
number of basis functions
is therefore
smaller, about 2,000.
And in MP2, you actually require
more basis functions per atom
to get an accurate
answer, so that the number
of feasible atoms is more
like 50 rather than 1,000.
So that's the cumulative effect
of larger basis sets and worse
scaling with number of atoms
that makes MP2 deal with much
smaller systems.
So questions about that.
Yes.
AUDIENCE: So what are some
examples of [INAUDIBLE]??
TROY VAN VOORHIS:
So C60 is 60 atoms.
That's an easy one.
There are a number
of small dyes,
for example, that we worked with
in our group, where you might
have seen absorption
spectra, or emission spectra,
or HOMOs, or LUMOs,
or hole transport
properties, or electron
transport properties,
that are around
50 atoms in size.
And then the main thing that
people get interested in that
are on the 1,000-atom regimes
are enzymes and peptides.
Those are the sort of--
when you start saying, I
want to do 1,000 atoms,
it's usually because it's really
some polymer or heteropolymer.
Another case where you also are
sometimes interested in things
where you do simulations
with 1,000 atoms,
if you're interested
in a surface.
Because if you have a
chunk of gold, for example,
1,000 atoms is just 10 gold
atoms, by 10 gold atoms,
by 10 gold atoms, which is just
a little chunk-- a very, very
small chunk of gold,
but still a chunk.
Other questions.
OK, so I'm going to spend seven
minutes talking about density
functional theory.
And then we're going to go
over and do some examples.
So the idea of density
functional theory
is that it'd be really nice
if you had some magical way
to do a Hartree-Fock
calculation,
but have a give you
exactly the right answer.
That basically would be the
dream because we saw, well,
for Hartree-Fock, we can
do big systems, it's cheap,
it has good scaling, all
these kinds of things.
It just doesn't give
us very good answers.
The results are pretty poor.
So what if we had
something that scaled
like it had the computational
cost of Hartree-Fock,
but was, say, the accuracy of
MP2, or even better than that?
And so in density functional
theory, what we use
is the idea of looking at
the electron density rho of r
as the fundamental variable.
So you can actually work
out, for a determinant,
what the electron density is.
And it's just the
sum of the squares
of the orbital densities--
or the sum of the
orbital densities.
So you square each orbital,
then you add that up,
and that gives you the
total electron density.
It's basically just
the probability
of finding an electron
at a point, r.
OK, so that's a nice observable.
The thing that's
kind of amazing is
that there is a theorem, a
mathematical theorem that you
can prove in about two pages
using just sort of proof
by contradiction.
You can prove that there exists
a functional, which is always
given then in e
sub v of rho, such
that when you're given the
ground state density, rho
0, if you plug that into
this magical functional,
it will give you the
ground state energy.
So if you found the ground
state density lying around
in the gutter, and
you picked it up,
and you put it into
this functional,
it would give you the
ground state energy.
And e0 is not just the
approximate ground state
energy.
It is the exact ground state
energy, the exact thing.
Further, for any
density, rho prime,
that is not the
ground state density,
if you plug that density into
eV, you get a higher energy.
So what you could do is
you could say, well, let me
just search.
Let me try density, see
what energy it gives.
Then I'll try another density.
If it gives a lower energy,
that's a better density.
And I'll keep going, and going,
and going, and going, until I
find that ground state density.
So the idea is that
then you would say, aha,
if I solve that equation--
so I search over densities,
and that's not going to
be a very hard search,
because the density is just
a three-dimensional object.
It's not like the wave
function that depends
on a bunch of coordinates.
It's just three dimensional.
So I solve that equation,
that will give me rho 0.
And then I can take that rho
0, plug it back in there,
and I will get the
ground state energy.
That gives me a closed
set of conditions
that lets me find the
ground state energy,
and then report back the
energy of that density.
AUDIENCE: Where do we
get this [? EV ?] from?
TROY VAN VOORHIS:
That's a great question.
The proof is a proof
by contradiction.
It's not constructive.
It proves that it
exists, but gives you
no way of constructing it.
So kind of the frustrating thing
is like, oh, cookies exist,
but we don't know
how to make them.
But the idea that
such a thing exists
gives you hope to
say, well, maybe
we can construct-- if we
can't find the exact one,
maybe we can find one
that's very, very good.
And that one that's very
good, we would just use.
We would pretend it was
exact, minimize the energy,
find the density, and then
report back the ground
state for that approximating.
And that's actually what you do
in density functional theory.
You have approximate
functionals.
And they're all-- just
like with basis sets--
named after the authors of the
people who wrote the papers.
Almost all of them are.
So there's the Local Spin
Density Approximation, LSDA.
There's the
Perdew-Burke-Ernzerhof
functional.
There's the
Becke-Lee-Yang-Parr functional,
the Perdew-Burke-Ernzerhof
Zero functional,
the Becke 3
Lee-Yang-Parr functional,
and then on and on and on.
Now these things
have been sort of--
there's a mishmash of
different exact conditions
that were used to derive
these functionals.
And then, empirically,
people have gone and shown
that PBE is better than LSDA.
It's about as good as BLYP.
BLYP is not as good as PBE0,
but it was just about as good
as B3LYP.
So we have, then,
sort of an idea
that, OK, if we go over towards
the right-hand side of this,
we're going to get better
results out of DFT.
But DFT, because it's based
off of a Slater determinant,
has exactly the same
computational scaling,
and storage, and
everything, as Hartree-Fock.
It's the same kind of expense.
But if I go back to
this accuracy thing,
and I say, well, what
would I get from that best
density functional
B3LYP, well what
I find is that activation
energies are still a little bit
difficult. Now I underpredict
activation energies.
But my dipole moments are good.
My sizes are again accurate.
My absorption energies are
just as accurate as MP2,
and my bond lengths are
just as accurate as MP2.
There are other things,
but basically B3LYP
is as accurate or
more accurate than MP2
for virtually every property
you would want, but has
Hartree-Fock-like cost.
And so that makes
density functional theory
the workhorse of computational
chemistry these days.
Because it's inexpensive, you
can do lots of calculations
with it, but it's accurate
enough for grunt work.
And then there's a blank column
there that, if I had time,
I would talk about.
There's a whole other
category of things
where you say, well,
perturbation theory doesn't
work, but are there
more constructive ways
that I could use wave functions
to approximate the correlation
energy?
And the answer is yes,
you can get good energies
and good answers out.
But then it comes at the cost
of the calculation getting
much more expensive.
So there are methods that
improve all these properties
on the wave function side,
that just require more.
So if you want to
know more about those,
they're in the notes.
And so now I think I will
switch over and show you how we
do some calculations with this.
So now I will switch.
There's two different
tools that you
have available to you
for free on Athena
that you can use to do
electronic structure
calculations.
One is Gaussian, and
the other one is Q-Chem.
There may be even
some other ones,
but those are the two
that I've collected notes
about how to use.
So I'm going to show
you how we use Q-Chem.
I'm not going to
do it on Athena.
I'm going to do it on my laptop,
but it's basically the same.
So there's this
GUI called IQmol.
And when you open
it up, it gives you
a window that looks like this.
And there's a few things here.
So there's, like, an open
file, you can create something,
you can open a file,
you can save something,
you can move something
around after you've built it.
And then these are the
build tools over here.
So that just turns
it into build mode.
If I click on this, I can choose
any item in the periodic table.
Right now, it's on carbon.
This lets me select
various fragments,
and also what the attach
point of the fragment's
going to be in some cases.
And I don't want to do it,
so I'll select it for now,
but then I'm going to
go back to this one.
This is the Add
Hydrogens button.
So if you have something that--
you just wrote your
Lewis structure,
it had no hydrogens in
it, you click that button,
it'll put all the hydrogens
where they should be,
or where it thinks
they should be.
This is a Minimize
Energy thing, which
basically if your structure
looks really crazy,
it'll kind of make
it look less crazy.
And then this is the
Select button, which
lets you pick certain atoms.
This lets you delete certain
atoms, take a picture,
do a movie.
And then that changes it
over to fullscreen mode.
And that's the life
preserver for Help.
But we'll do
building a molecule.
So does somebody want to
tell me a molecule that
has less than 15
atoms in it that they
would like me to draw here,
that I can actually draw?
AUDIENCE: Diethyl ether.
TROY VAN VOORHIS: Diethyl ether.
Yes, I can do that.
OK.
All right, so put
a carbon, carbon,
then I got to switch to an
oxygen, oxygen, carbon, carbon.
And then I will click
the Add Hydrogen button.
[LAUGHTER]
There we go.
I've got my diethyl ether now.
And while my geometry
here is not perfect,
it doesn't look totally crazy.
So I'll just go with that.
And then the thing
we're going to be
most interested in is
this little tab here,
called Calculation.
So calculation here has a little
button that says Q-Chem Setup.
So I'll do Q-Chem Setup.
And now it's got a bunch of
things that I can specify.
There are a few things
that are most important.
So the first thing is, what kind
of calculation I want it to do.
So I could have it
just compute the energy
of the molecule for the
configuration I specified.
That would be one useful thing.
I could compute forces.
In other words, the
forces on the atoms,
where the atoms want to move.
I could optimize the
geometry, which would then
relax the thing down, and figure
out what the best geometry is.
I could do various
things that scan
across the potential
energy surface.
I can compute
vibrational frequencies.
I can compute NMR
chemical shifts.
I can do molecular dynamics.
I can compute properties--
so lots of things.
I'm just going to
optimize the geometry.
That's going to be what I'm
going to choose for now.
And then I get to
choose a method.
So I can do Hartree-Fock.
I can do B3LYP.
I can go down-- so all the
things above this dotted line
here are different
density functions.
And then, down here, we've
got MP2, MP2 bracket v--
I'm not actually
sure what that is.
I've got various versions of
MP2 that try to make it faster.
And then I've got a
couple of cluster methods
that we didn't talk about yet,
but that are in the notes.
And then various
methods for getting
electronic excited states,
down here at the bottom.
So for now, let's just
do B3LYP, and then I
get to choose a basis.
Because this is on my
laptop, and we're in class,
I'm going to choose
a small basis.
So I'm going to leave it at
321g, which is a small basis.
And then I can also
specify the charge, which
is like the total
charge of the molecule,
and the spin multiplicity.
So I want it to be a
neutral molecule, spin 1--
or spin multiplicity
1, which is spin 0.
Because it's 2s plus 1
is the spin multiplicity.
And then there's various things
about convergence control
that you probably won't
need to worry about.
And then there's advanced
things that you also probably
won't need to worry about.
So now I'm happy with that.
I can submit my job.
And I will call this
diethyl ether, hit OK.
And now it doesn't appear
to be doing anything.
In the background, it
is doing something.
I think this will finish
in just a few seconds,
but if you get concerned that
the job might not be running
or something might
have happened,
you can go over to
the Job Monitor.
Methane-- that was the job
that I ran this morning
to make sure that it works.
Diethyl ether, we just submitted
it at 10:50 and 21 seconds.
Ah, there we go.
Now it's been running
for four seconds.
Let's see how long it takes.
OK, it's taking a lot
longer than I think.
All right.
AUDIENCE: They close
[? in 15 ?] [INAUDIBLE]..
TROY VAN VOORHIS:
That's all right.
We still we've
got a few minutes.
So what will happen here
in a minute, when we--
AUDIENCE: [INAUDIBLE].
TROY VAN VOORHIS:
That's all right.
It's going along.
So when it's done,
what will happen
is, over here,
it'll let us know.
And over here will
appear a thing
where we can start
looking at the results.
It might have gone faster if
I had minimized the energy
to begin with, because then it
would do less geometry steps.
It's going along.
I don't know.
We might run out of class
time before it finishes.
It's 10:51.
We've got four minutes.
Let's see.
Come on, laptop.
I've got the MacBook Pro.
Maybe I should have upgraded the
processor, this kind of thing.
AUDIENCE: How long
do you usually
let your group
calculations go for?
TROY VAN VOORHIS: Well,
in my research group,
we have 1,000 cores that we--
maybe around 1,200 cores,
even, I think, actually.
Oh, there we go.
It finished.
But we have around 1,200
cores for around 10 people.
So that means that you can let
a job run for a very, very, very
long time, on
multiple cores, and it
doesn't get in anyone's way.
So it's not unusual for us to
put something in there that
lasts a week.
So now it's finished.
We want to copy the
results from the server.
This is an arcane
way of saying things.
Because it's designed
for supercomputers, where
you are sitting at a
terminal in your office,
and then it sends
all the results over
to the supercomputer,
it runs the calculation,
and then you have to
copy the results back
to your little computer.
I was actually running
it on my computer.
It's just the computer is both--
I'll just put it--
New Folder.
OK.
There we go.
Now it has a little star
next to it, which lets
me know that it's finished.
So let's see, Info--
so I can show the dipole moment.
So there we go, that
little blue arrow--
which, magnitude and
direction is the dipole moment
of diethyl ether.
It does have a dipole moment,
with its positive end here,
pointing in that direction.
That seems about right.
Let's see.
I can look at the files.
These are only if you're really,
really interested in details.
You see these are very verbose.
They have every
bit of information
you could possibly want in them.
So if there's ever something
that you can't get out
of the GUI, if you
just kind of go--
you can search.
And I can look for converge.
Oh, look, optimization
converged, and then it
tells me everything after that.
I can look at the different
atoms in the molecule,
and it will highlight
them for me.
And then I can look at
bonds, and it will tell me--
I don't if you guys can
read, but down here, it
has the bond length
in the corner.
So that C-H bond is 1.093.
All these C-H bonds are going
to be 1.09, plus or minus a bit.
The carbon-oxygen bond
there is 1.46, rounding.
That one is also 1.46.
And then this here-- so as it
did the geometry optimization,
it was adjusting the geometry
trying to find the minimum.
And these are the
different total energies
that it had as it went
through the optimization.
You'll see that eventually
it sort of slows down.
The energy is going
down every time,
and eventually it sort
of slows down to where
it's not changing any more.
It says, all right,
that's the minimum.
You also notice that
these are in atomic units.
So that's minus 232 Hartrees.
So that's like 4,000 kcals.
No, that's more
than 4,000 kcals.
Anyway, it's tens of
thousands of kcals.
The reason it's
such a low number
is this is the energy it would
take to rip every electron off
of the molecule, including
all of the 1s electrons.
So the energies
are typically going
to be very, very large
negative numbers.
Don't worry about that.
It's energy differences
that matter.
So you notice, even my bad
geometry was already 232.1.
The correct answer is 232.3.
So the total energy
change wasn't that much.
And then the last
cool thing that I
want to show you is that
we can plot orbitals.
So I'm going to
count some orbitals.
There we go.
So I just had to calculate
the HOMO and the LUMO.
I didn't go through
all of it with you.
But if you go there, it
gives you some things--
I should actually go
back and look, show you.
So you get to choose
a range of orbitals.
The iso value is basically
where it draws the surface.
So it's inverse of
what you expect.
A larger iso value is--
it's an isosurface,
so if you make it larger,
the place where the iso
value is larger is closer in.
So it makes a smaller bubble.
And a smaller iso value
makes a bigger bubble.
And you can change
the number of points
that it samples to make the
surface, change the colors,
change the opacity
of it in terms
of how you want it to look.
But I already did this.
So there is the HOMO.
So the HOMO is something
in the pi manifold.
It's pi star here.
And then I can take out
the HOMO and show the LUMO.
And the LUMO is, in this
particular basis set,
not very useful, because I don't
have a lot of basis functions.
So my basis functions, just
there's not enough of them
to really describe the
LUMO, or the LUMO plus 1.
If I want to get
a LUMO, we'd need
a bigger basis set than this.
So that's the kind of
thing that you can do
with Q-Chem or with Gaussian.
Gaussian has a GUI
called GaussView.
You can use either one.
But they both work pretty
well, and are pretty intuitive.
And so I think, on
the next problem set,
you will have some
homework problems
that involve running some
calculations with these things.
And I hope you enjoy it.
Well, I don't know if homework
is ever really enjoyable.
But I hope it's marginally
enjoyable anyway.
