Implementation of pair potential methods

Implementation of pair potential methods

(Parte 1 de 2)

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 1

Material for these lecture notes was compiled from the references below

9MIT’s 3.320 course notes (Prof. G. Ceder)

9A.P. Leach, Molecular Modeling: Principles and Applications

9Daw, M. S., Foiles, S. M. & Baskes, M. I. The EAM: a review of theory and applications Materials Science Reports 9, 251 (1993).

9Example of the Clementi and Roetti Tables 9Collection of papers on empirical potentials for Silicon

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 2

Pair potential methods

We will discussto some detail about the failures of pair potential modelsand how to address them with other empirical models.

Recall that in a pair potential model, your energy is written pairwise, i.e. for each pair of atoms (i,j) you need to evaluate the distance between them and, then, use that as an argument in some kind of energy function.

This is an N2 operationwhere Nis the number of atoms.

You want to use pair potentials(e.g. rather than QM-methods) when your system gets big. But you need to do better than N2 .

If you want to calculate energies or forces on an atom you only need to account for interactions with atoms within a certain distance. Thus you keep neighbor lists, e.g. for every atom you figure out the atoms that are possibly in its neighbor environment.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 3

Pair potential methods

For solids, coming up with neighbor lists is easy--the topological relation between atoms would not change much during the simulation.

For a gas or a liquid you need to regularly update the neighbor lists.

Even though there is an overhead with this update, the overall cost can be reduced to order Nrather than N2 .

We usually want to compute the lowest energy of the system (relaxation). To do so (using MD), we need to compute the forces to figure out how atoms accelerate or decelerate.

Computing the forcesrequires computing the derivative of the total energy with respect to a given position –that is an N2 operation(for each of the N-atoms, you have to sum over pairs of atoms).

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 4

Pair potential methods

The minimization is done by conjugate gradient methodswith line minimization, etc. (see book by A.P. Leach for more details).

With known force and assuming you know what the curvature is of the energy in that direction you can move towards its minimum.

This is an iterative process. It has been shown to work very well for very large systems.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 5

There is a notch on the top and bottom. You pull from the side and you see the crack growing and dislocations sprout (they counted under-or over-coordinated atoms). With 1 billion atoms things start looking like a continuum.

Large scale atomistic simulation

A million atoms isa cube of size ~ 200 Angstroms!So these are still very small systems!

IBM Almaden

This is a billion atom MD simulation. Dislocation lines are shown.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 6

Periodic boundary conditions

To get rid of boundary effects you use periodic boundary conditions in simulations and in as many directions as you can.

With periodic boundary conditions if you have some unit (Fig. above) you, essentially repeat it next to it. In principle you manage to have infinite atoms in your simulation but finite degrees of freedom!

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 7

Periodic boundary conditions

How do we figure out what repeat unit to use?

¾For a periodic material (e.g. a perfect crystalline material) like the peroxide unit cell shown there is no problem or approximations to be made (at least for static calculations).

¾When you do dynamics, the translational symmetry can be temporarily broken because each atom can vibrate independent of other atoms!

In the dynamic case, imposing periodicity implies that you assume that all atoms vibrate together(you are imposing a cut-off on the frequencies of vibrations that can occur).

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 8

Periodic boundary conditions and defects

¾If you take a periodic structure and you introduce a defect (Oxygen vacancy), periodicity is broken.

In this case, you can impose periodicity by repeating the defect

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 9

Supercell approximation

You now have a bigger unit cell with the defect in it, here an Oxygen vacancy (supercell approximation).

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 10

Supercell approximation

You need to make a cell that is big enough so that the defects do not interact, because, then, if you have Ncells, your calculation just has the energy of N defects.

If the defects interact, then, you are not modeling an isolated defect.

You should always check convergence related to the size of the supercell(e.g. use bigger and bigger supercelland if the answer does not change you are done).

Rarely you have potentials that are relevant over large distances that you need to make really large supercells.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 1

Supercell approximation –Long range effects

However, you may have other energetic effects (long range) ¾electrostatics(that gives you very long range interaction)

¾charge defects:If we did a charge defect in the supercell of peroxide (e.g. we took out an Oand the supercell has now a net positive charge), you will have an infinite energy if you do an infinite system(positive charge interacting with only other net positive charges).

(the 2 electrons ofOwill go some-where), then you will

If you decide to compensate with a negative charge have dipoles interacting.

Dipoles will not give you infinite energy but they will give you a fairly long range energy effect.

You can often do much smaller cell calculation by explicitly taking out the defect-defect interaction.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 12

Supercell approximation -relaxation

Often causing problems isrelaxation: e.g. if you put a big ion in the center e.g. of the peroxide unit cell, you would build up a strain field. The defects would interact through their strain fields.

Another way is when the volume of a supercellas it relaxes interferes with the relaxation of the supercellnext to it (homogeneous strain interaction)

In conclusion, you rarely interact through the direct energetics, you often interact through secondary effects

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 13

Calculation of the vacancy formation energy in aluminum

The figure shows the vacancy formation energy versus the number of atoms in the supercell(quantum mechanical simulation)

Note that from this figure is quite hard to figure out whether you have converged.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 14

Failures of pair potentials

It became clear in the 80’s that there are certain things you could not get right with pair potentials because of problems of the wrong form’ rather than of wrong parameters’.

Typical example of such failure was the vacancy formation energy calculation.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 15

Failures of pair potentials: Vacancy formation

Let us calculate the pair approximation to the vacancy vacancy formation energy for FCC crystal.

How do we make a vacancy?You take the central atom and you move it somewhere else in the bulk

In making the vacancy, we destroy 12 bonds because this atom we remove is coordinated 12-fold in FCC (thus we loose 12 times the bond energy).

Of course we put this energy back somewhere in the bulk and so we gain the cohesive energy of one atom.

What’s the cohesive energy per atom in Fcc?It’s 6 times the bond energy because there are 12 bonds around an atom but every bond is shared between 2 atoms. So, there are 6 per atom.

ΔEvacancy= -12 E bondenergy+ 6 E bond energy =

= -6 E bond energy = -E cohesive energy /atom

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 16

Vacancy energies So this is what you find with pair potentials:

This is a static approximation–we assumed that all bonds were unchanged when we took the atom out.

In reality, if I take that central atom out, the atoms around it will relax somewhat and change their bond energy.

ΔE vacancy= -12 Ebond energy+ 6 Ebond energy =

= -6 Ebond energy

= -Ecohesive energy /atom

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 17

Cohesive energy

Let us look at some results of the cohesive energy. For noble gases you’re reasonably close

Note that noble gases are just inert shells interacting through van der Waals interactions and that’s how we came up with the LJ potential.

For metals, the vacancy formation energy is, actually, only a small fraction of the cohesive energy (can never get that low with pair potentials) –the reason is because that there is certain physics missing from pair potentials.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 18

Cohesive energy

Why is the vacancy formation energy so low in real materials?

Remember that in metals the bonding energy does not go linear in the coordination number (Z) but like the square root of Z.

When you take out an atom, you make a lot of the atoms around it go from

12 coordinationsto 1 and they loose 1/12th of the cohesive energy because the cohesive energy is the average of all bonds.

Let us make a simple model by taking the cohesive energy as

E cohesive energy

If we put the atom somewhere back in we would gain Cand, then, as

Let us calculate the vacancy formation energy in this model, for FCC (Z =12).

we take it out, what do we lose? We go from coordination 12 to 1. 12 vacancy co h

So, in a square root model, you’re already closer to the experiments where the vacancy formation energy is a much smaller fraction of the cohesive energy .

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 19

Surface relaxation

If you create a surface in a metal, the distance between the 1st and the

2nd layers becomes smaller. So, in a metal, the 1st layer tends to relax inward.

In pair potentials, the 1 stand the 2nd neighbor distance is often bulk like and, in some cases, actually, it relaxes outwards because of competing forces. But, in real systems, it tends to relax inward. Why is that?

The surface plane relaxes inwards because the bonds between the 1st and the

2nd layers strengthen because the surface is under-coordinated.Although you are loosing energy (that’s why there is surface energy), the energy per bond increases.

Remember, if you go to lower coordination, the energy per bond is, actually, becoming stronger. If you cut away the bonds (blue) to make the surface, the remaining bonds (red) strengthen and they pull the layer inwards.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 20

The Cauchy problem

The Cauchy problem (relation between the stress and the strain tensor) is another problem that is common in pair potentials.

For pair potentials C12and C44 are always the same. This has to do with shearing being a zero energy strain mode (see last lecture).

In a cubic model you can slide it without any energy change and this

For metals: C12/C44 =1.5 for Cu, 1.9 for Ag, 3.7 for Au. For van derWaals solids and ionic crystals, the Cauchy relation is satisfied C12/C44 =1.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 21

Predicting crystal structures with pair potentials

You cannot predict crystal structures in any thing that is covalent (metals are delocalized covalent).

Pair potentials tend to give you close packingbecause they want to have as many atoms around as possible in order to bring the cohesive energy down.

Also pair potentials do not count that if you bring more atoms around your bond strength weakens.

Thus they miss these 2 important components to get crystal structure right.

You need at least the fourth moment of ``the density of states’’ to get crystal structure right (a four-body effect).

This means that if we write the energy as a constant plus an expansion in two-body term, and an expansion in a three-body term, etc. (cluster expansion), we will need to include at least a four-body term to converge.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 2

How to fix pair potential problems

You can make cluster potentials(used a lots in organic materials, polymers, etc.)

¾If you do not like a two-body potential, you can make e.g. a three-body potential (this will account for angular dependence, bonding angles)

An alternative is to use pair functionals(used in metals). Here, you still do pairwise summation but the cohesive energy per atom does not go linear with what you sum pairwise.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 23

There is a whole class of methods, which are called effective medium theories

Effective Medium Theories The embedded atom method (EAM) (Mike Baskes) is the most popular.

You start from the understanding that the energy should not be linear with coordination.

You write the energy per atom as a function of the number of bonds around it

The function fshould not be linear.

These are called energy functionalsrather than pair potentials. How do you measure the number of bonds?

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 24

EAM: Measuring the number of bonds

¾You dontwant to do this discretely. The problem is, as soon as you measure bonds discretely within some radius, you’re going to have to deal with discontinuities. As atoms move in and out of that radius, you got large forces because of these discontinuities around your cut-off interface.

¾You need a smooth count of coordination. In EAM this is based on measurement of the electron density being projected on you from your neighbors(the book on electronic structure of materials by A.P. Sutton does a good presentation on this).

As atoms move further away, they still contribute but they contribute less. So, this is a continuous measure.

The more electron density you project, either the closer-by the atoms are or the more of them there are.

How do you measure the number of bonds?

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 25

Atomic electron densities

You calculate the electron densityon a site as a measurement, in some sense, of the coordination, as a sum of electron densities coming from its neighboring atoms j.

atomic electron density of atom j

In EAM, atomic electron densities used are usually computed with Hartree-Fockmethods. They are parameterized by Gaussians.

Atomic densities are tabulated in E. Clementiand C. Roetti , Atomic Data and

Nuclear Data Tables, Vol14, p. 177 (1974).

You treat the electron density of the crystal as the overlap of the atomic densities of the atoms

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (3/30/2009) 26

Energy in EAM

There are two parts in the energy in EAM. At first, you have the embedding energy. To compute it, you sum over every atom and you are looking for its cohesive energy, which is some function Fof the projected densities of all the other atoms j, which is the electron density at the site of i.

This is the non-linear part of the energy.

We have to solve 3 problems: ¾The electron densities are taken from the Clementi & Roetti tables.

¾We need to know what the embedding function is.

¾The pair potential used can be almost anything (it’s mainly used so the atoms don’t fly into each other). Typically, they use some screened electrostatic form , with . Even LJ will be OK.

(Parte 1 de 2)