Loading in 2 Seconds...

Ch121a Atomic Level Simulations of Materials and Molecules

Loading in 2 Seconds...

- 87 Views
- Uploaded on

Download Presentation
## PowerPoint Slideshow about 'Ch121a Atomic Level Simulations of Materials and Molecules' - yardan

**An Image/Link below is provided (as is) to download presentation**

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript

### Ch121a Atomic Level Simulations of Materials and Molecules

### Ch121a Atomic Level Simulations of Materials and Molecules

Lecture: Monday, Friday 2-3pm

Lab Session: Wednesday 2:30-3:30pm

Extra help: Monday 3-3:30pm and Wednesday 2-2:30pm

Lecture 3, April 8, 2013

FF1: valence, QEq, vdw

William A. Goddard III, wag@wag.caltech.edu

Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology

Hours: Monday, Wednesday, Friday 2-3pm

Lecture 3, April 8, 2013

FF1: valence, QEq, vdw

William A. Goddard III, wag@wag.caltech.edu

Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology

TA’s Caitlin Scott and Andrea Kirkpatrick

Homework and Research Project

First 5 weeks: The homework each week uses generally available computer software implementing the basic methods on applications aimed at exposing the students to understanding how to use atomistic simulations to solve problems.

Each calculation requires making decisions on the specific approaches and parameters relevant and how to analyze the results.

Midterm: each student submits proposal for a project using the methods of Ch121a to solve a research problem that can be completed in the final 5 weeks.

The homework for the last 5weeks is to turn in a one page report on progress with the project

The final is a research report describing the calculations and conclusions

The full Quantum Mechanics Hamiltonian: Htotal(1..Nel,1..Mnuc)

Hel(1..Nel) =

N

M

M,N

M

N

nuclei

electron-nucleus

electrons

nucleus-nucleus

electron-electron

Born-Oppenheimer approximation, fix nuclei (Rk=1..M), solve for Ψel(1..N)

Hel(1..Nel) Ψel(1..Nel)= EelΨel(1..Nel)

Ψel(1..Nel)and Eel(1..Mnuc) are functions of the nuclear coordinates.

Next solve the nuclear QM problem

Hnuc(1..Mnuc) Ψnuc(1..Mnuc)= EnucΨnuc(1..Mnuc)

Hnuc(1..Mnuc) =

+Eel(1..Mnuc)

FF is analytic fit to this

EPE(1..Mnuc)

The Bending potential surface for CH2

C

H

H

1B1

1Dg

1A1

3B1

3Sg-

9.3 kcal/mol

linear

Note that E(p-) = E(p+)

Symmetric about = p=180°

EPE(1..Mnuc)

= Eel(1..Mnuc)

The Force Field (FF) is an analytic fit to such expressions but formulated to be transferable so can use the same parameters for various molecules and solids

General form:

Ecross-terms

Involves covalent bonds and the coupling between them (2-body, 3-body, 4-body, cross terms)

E nonbond = Evdw + Eelectrostatic

Involves action at a distance, long range interactions

Bond Stretch Terms: Harmonic oscillator form

R

Taylor series expansion about the equilibrium, Re, d = R-Re

E(d)=E(Re)+(dE/dR)e[d]+½(d2E/dR2)e[d]2+(1/6)(d3E/dR3)e[d]3 + O(d4)

ignore

ignore

zero

ignore

Kb

Re = 0.9 – 2.2 Å

Ke = 700 (Kcal/mol)/Å 2

E(R) = (ke/2)(R-Re)2,

Simple Harmonic

Units: multiply by 143.88 to convert mdynes/ Å to (Kcal/mol)/Å 2

QM eigenstates, En = n + ½ , where n=0,1,2,…

Ground state wavefunction (Gaussian)

0(d) = (mw/pЋ) exp[- (mw/2Ћ)d2]

Re

Some force fields, CHARRM, Amber use k=2Ke to avoid multiplying by 2

Problem: cannot break the bond, E ∞

Bond Stretch Terms: Morse oscillator form

R

Re

De

Want the E to go to a constant as break bond, popular form is Morse function

where X = exp[-a(R-Re)] = exp[-(g/2)(R/Re -1)]

Evdw (Rij) = De[X2-2X]

Note that E(Re) = 0, the usual convention for bond terms

De is the bond energy in kcal/mol

Re is the equilibrium bond distance

a is the Morse scaling parameter

calculate from Ke and De

- = sqrt(ke/2De),

where ke = d2E/dR2 at R=Re

(the force constant)

E = 0

Morse Potential – get analytic solution for quantum vibrational levels

where X = exp[-a(R-Re)] = exp[-(g/2)(R/Re -1)]

EMorse(Rij) = De[X2-2X]

a = sqrt(ke/2De), where ke = d2E/dR2 at R=Re (the force constant)

Ev = hv0(n + 1/2 ) – [hv0(n + ½)]2/4De

n0 = (a/2p)sqrt(2De/m)

Level separations decrease linearly with level

En+1 – En = hv0 – (n+1)(hv0)2/2De

Write

En/hc = we(n + ½ ) - xewe(n + ½ )2

D0

De

Angle Bend Terms: cosine harmonic

J

K

I

Expand E() as Fourier series.

Note only cos(n) since must

be symmetric about =0 and p

E() = C0 + C1cos() + C2cos(2) + C3cos(3) + ….

E() ~ ½ Ch [cos - cos e]2

E=0 at = e

E’() = dE()/d = - Ch [cos - cos e] sin = 0 at = e

E”() = d2E()/d2 = - Ch [cos - cos e] cos + Ch (sin )2

k = E”(e) = Ch(sin e)2

Using cos(2) = 2 cos2 - 1

The Bending potential surface for CH2

C

H

H

1B1

1Dg

1A1

3B1

3Sg-

9.3 kcal/mol

linear

Note that E(p-) = E(p+)

Symmetric about = p=180°

Barriers for angle term

The energy is a maximum at = p.

Thus energy barrier to “linearize” the molecule becomes

By symmetry the angular energy satisfies

This is always satisfied for the cosine expansion but

A second more popular form is the Harmonic theta expansion

E() = ½ K [ - e]2

However except for linear molecules this does NOT satisfy the symmetry relations, leading to undefined energies and forces for = 180° and 0°. This is used by CHARMM, Amber, ..

Angle Bend Terms for linear molecule,

K

J

I

If e = 180° then we write E() = K [1 + cos()]

since for = 180 – d this becomes

E(d) = K [1 + (-1 + ½ d2)] ~ 1/2 Kd2

Angle Bend TermsSimple Periodic Angle Bend Terms

For an Octahedral complex

The angle function should have the form

E() = C [1 - cos4] which has minima for

= 0, 90, 180, and 270°

With a barrier of C for

= 45, 135, 225, and 315°

Here the force constant is K= 16C = CP2where P=4 is the periodicity, so we can write C=K/P2

Angle Bend TermsSimple Periodic Angle Bend Terms

- For an Octahedral complex
- The angle function should have the form
- E() = C [1 - cos4] which has minima for

= 0, 90, 180, and 270°

With a barrier of C for

= 45, 135, 225, and 315°

Here the force constant is K = 16C = CP2 where P=4 is the periodicity, so we can write C=K/P2

However if we wanted the minima to be at

- = 45, 135, 225, and 315° with maxima at

= 0, 90, 180, and 270°

then we want to use E() = ½ C [1 + cos4]

Thus the general form is

E() = (K/P2)[1 – (-1)Bcos4]

Where B=1 for the case with a minimum at 0° and B=-1 for a maximum at 0°

Dihedral barriers

7.6 kcal/mol Cis barrier

HOOH

Part of such barriers can be explained as due to vdw and electrostatic interactions between the H’s.

But part of it arises from covalent terms (the pp lone pairs on each S or O)

This part has the form E(φ) = ½ B [1+cos(2φ)]

Which is

E=0 for φ=90,270° and E=B for φ=0,180°

φe ~ 111°

1.19 kcal/mol Trans barrier

HSSH: φe ~ 92°

5.02 kcal/mol trans barrier

7.54 kcal/mol cis barrier

dihedral or torsion terms

J

L

φ

L

I

J

K

K

I

Given any two bonds IJ and KL attached to a common bond JK, the dihedral angle φ is the angle of the JKL plane from the IJK plane, with cis being 0

E(φ) = C0 + C1 cos(φ) + C2 cos(2 φ) + C3 cos(3 φ) + ….which we write as

Where each kn is in kcal/mol, n = 1, is the periodicity of the potential and

d=+1 the cis conformation is a minimum

d=-1 the cis conformation is thea maximum

Input data is kn and dnfor each n.

Consider the central CC bond in Butane

L

I

J

K

There are nine possible dihedrals: 4 HCCH,

2 HCCC, 2 CCCH, and 1 CCCC.

Each of these 9 could be written as

E(φ) = ½ B [1 + cos(3 φ)],

For which E=0 for φ = 60, 180 and 300°

and E=B the rotation barrier for φ = 0, 120 and 240°.

Consider the central CC bond in Butane

L

I

J

K

There are nine possible dihedrals: 4 HCCH,

2 HCCC, 2 CCCH, and 1 CCCC.

Each of these 9 could be written as

E(φ) = ½ B [1 – cos(3 φ)],

For which E=0 for φ = 60, 180 and 300°

and E=B the rotation barrier for φ = 0, 120 and 240°.

For ethane get 9 HCCH terms.

The total barrier in ethane is 3 kcal/mol but standard vdw and charges of +0.15 on each H will account for ~1 kcal/mol.

Thus only need explicit dihedral barrier of 2 kcal/mol.

Consider the central CC bond in Butane

L

I

J

K

There are nine possible dihedrals: 4 HCCH,

2 HCCC, 2 CCCH, and 1 CCCC.

Each of these 9 could be written as

E(φ) = ½ B [1 + cos(3 φ)],

For which E=0 for φ = 60, 180 and 300°

and E=B the rotation barrier for φ = 0, 120 and 240°.

For ethane get 9 HCCH terms.

The total barrier in ethane is 3kcal/mol but standard vdw and charges of +0.15 on each H will account for ~1 kcal/mol.

Thus only need explicit dihedral barrier of 2 kcal/mol.

Can do two ways.

Incremental: treat each HCCH as having a barrier of 2/9 and add the terms from each of the 9 to get a total of 2 (Amber, CHARRM)

Or use a barrier of 2 kcal/mol for each of the 9, but normalize by the total number of 9 to get a net of 2 (Dreiding)

Twisting potential surface for ethene

The ground state (N) of ethene prefers φ=0º to obtain the highest overlap of the pp orbitals on each C

The rotational barrier is 65 kcal/mol.

We write this as E(φ) = ½ B [1-cos2 φ)]. Since there are 4 HCCH terms, we calculate each using the full B=65, but divide by 4.

T

The triplet excited state (T) prefers φ =90º to obtain the lowest overlap.

N

φ = 90°

φ = 0°

φ = -90°

Inversions

L

J

I

K

L

J

I

K

When an atom I has exactly 3 distinct bonds IJ, IK, and IL, it is often necessary to include an exlicit term in the force field to adjust the energy for “planarizing” the center atom I.

Where the force constant in kcal/mol is

Umbrella inversion:

ω is the angle between

The IL axis and the

JIK plane

AMBER Improper Torsion JILK

J

L

I

K

Amber describes inversion using an improper torsion.

n=2 for planar angles (0=180°) and n=3 for the tetrahedral

Angles (0 =120°).

Improper Torsion: is the angle between the JIL plane and the KIL plane

CHARMM Improper Torsion IJKL

L

J

I

K

In the CHARMM force field, inversion is defined as if it were a torsion.

For a tetrahedral carbon atom with equal bonds this angles is φ=35.264°.

Improper Torsion: φ is the angle between the IJK plane and the

LJK plane

Summary: Valence Force Field Terms

example

Description

Typical Expressions

Harmonic Stretch

Bond Stretch

Harmonic-cosine

Angle bend

dihedral

Torsion

E(φ) = ½ B [1 – cos(3 φ)],

Umbrella inversion

Inversion

Coulomb (Electrostatic) Interactions

Most force fields use fixed partial charges on each nucleus, leading to electrostatic or Coulomb interactions between each pair of charged particles. The electrostatic energy between point charges Qi and Qk is described by Coulomb's Law as

EQ(Rik) = C0 Qi Qk /(e Rik)

where Qi andQk are atomic partial charges in electron units

C0 converts units:

if E is in eV and R in A, then C0 = 14.403

If E is in kcal/mol and R in A, then C0 = 332.0637

e = 1 in a vacuum (dielectric constant)

TA check numbers

How estimate charges?

Even for a material as ionic as NaCl diatomic, the dipole moment a net charge of +0.8 e on the Na and -0.8 e on the Cl not the idealized charges of +1 and -1

We need a method to estimate such charges in order to calculate properties of materials.

A side note about units.

In QM calculations the unit of charge is the magnitude of the charge on an electron and the unit of length is the bohr (a0)

Thus QM calculations of dipole moment are in units of ea0 which we refer to as au. However the international standard for quoting dipole moment is the Debye = 10-10esu A

Where m(D) = 2.5418 m(au)

QEq Electrostatics energy

I

atomic

interactions

Keeping:

Jij

Hardness (IP-EA)

Electronegativity (IP+EA)/2

1/rij

rij

ri0 +rj0

- Self-consistent Charge Equilibration (QEq)
- Describe charges as distributed (Gaussian)
- Thus charges on adjacent atoms shielded
- (interactions constant as R0) and include interactions over ALL atoms, even if bonded (no exclusions)
- Allow charge transfer (QEq method)

2

Three universal parameters for each element:

1991: used experimental IP, EA, Ri;

ReaxFF get all parameters from fitting QM

Charge dependence of the energy (eV) of an atom

E=12.967

E=0

E=-3.615

Cl+

Cl

Cl-

Q=+1

Q=0

Q=-1

Harmonic fit

The 2nd order expansion in Q is ok

Get minimum at Q=-0.887

Emin = -3.676

= 8.291

= 9.352

QEq: the optimum charges lead to Equilibrationthat is, equal chemical potential

Expand the energy of an atom in a power series of the net charge on the atom, E(Q)

Including terms through 2nd order leads to

- Charge Equilibration for Molecular Dynamics Simulations;
- K. Rappé and W. A. Goddard III; J. Phys. Chem. 95, 3358 (1991)

(2)

(3)

Interpretation of J, the hardness

Define an atomic radius as

RA0

Re(A2)

Bond distance of homonuclear diatomic

H 0.84 0.74

C 1.42 1.23

N 1.22 1.10

O 1.08 1.21

Si 2.20 2.35

S 1.60 1.63

Li 3.01 3.08

Thus J is related to the coulomb energy of a charge the size of the atom

The total energy of a molecular complex

Consider now a distribution of charges over the atoms of a complex: QA, QB, etc

Letting JAB(R) = the Coulomb potential of unit charges on the atoms, we can write

Taking the derivative with respect to charge leads to the chemical potential, which is a function of the charges

or

The definition of equilibrium is for all chemical potentials to be equal. This leads to

The QEq equations

Adding to the N-1 conditions

The condition that the total charged is fixed (say at 0) leads to the condition

Leads to a set of N linear equations for the N variables QA.

AQ=X, where the NxN matrix A and the N dimensional vector A are known. This is solved for the N unknowns, Q.

We place some conditions on this. The harmonic fit of charge to the energy of an atom is assumed to be valid only for filling the valence shell.

Thus we restrict Q(Cl) to lie between +7 and -1 and

Q(C) to be between +4 and -4

Similarly Q(H) is between +1 and -1

The QEq Coulomb potential law

We need now to choose a form for JAB(R)

A plausible form is JAB(R) = 14.4/R, which is valid when the charge distributions for atom A and B do not overlap

Clearly this form as the problem that JAB(R) ∞ as R 0

In fact the overlap of the orbitals leads to shielding

shielded

unshielded

The plot shows the shielding for C atoms using various Slater orbitals

Using RC=0.759a0

And l = 0.5

QEq for Ala-His-Ala

Amber charges in parentheses

QEq for deoxy adenosine

Amber charges in parentheses

Polarizable QEq

Allow each atom to have two charges:

A fixed core charge(+4 for Si)with a Gaussian shape

A variable shell charge with a Gaussian shape but subject to displacement and charge transfer

Electrostatic interactions between all charges, including the core and shell on same atom

Allow Shell to move with respect to core, to describe atomic polarizability

Self-consistent charge equilibration(QEq)

Four universal parameters for each element:

Get from QM

Noble gas dimers

All nonbonded atoms and molecules exhibit a very repulsive interaction at short distances due to overlap of electron pairs (Pauli Repulsion) and a weak attractive interaction scaling like 1/R6 at long R (London Dispersion. Together these are called van der Waals (vdW) interactions

Ar2

s

Re

De

Most popular form: Lennard-Jones 12-6

E=A/R12 –B/R6

= De[r-12 – 2r-6] where r= R/Re

Note that E is a min at Re, note the 2

- A 2nd form for LJ 12-6 is
- = 4 De[t-12 – t-6] where t=R/s

Note that E=0 at R = s, note the 1

Here s = Re(1/2)1/6 =0.89 Re

van der Waals interaction

- What do vdW interactions account for?
- Short range: Pauli Repulsion of overlapping electron pairs

Pauli repulsion:

Born repulsion:

- Long range attraction London dispersion due to Instantaneous fluctuations in the QM charges

Dipole-dipole

Dipole-quadrupole

Quadruole-quadrupole

We usually neglect higher order (>R6) terms.

Popular vdW Nonbond Terms: LJ12-6

NonBond

R

Lennard-Jones 12-6

E(R)=A/R12 – B/R6

=Dvdw{1/r12 – 2/r6} where r = R/Rvdw

Here Dvdw=well depth,

ndRvdw= Equilibrium distance for vdwdimer

Alternative form

E(R) = 4Dvdw{[svdw/R)12]- [svdw/R)6]}

Where s = point at which E=0 (inner wall, s ~ 0.89 Rvdw)

Dimensionless force constant k = {[d2E/dr2]r=1}/Dvdw = 72

The choice of 1/R6 is due to London Dispersion (vdw Attraction)

However there is no special reason for the 1/R12 short range form (other that it saved computation time for 1950’s computers)

LJ 9-6 would lead to a more accurate inner wall.

Popular vdW Nonbond Terms- Exponential-6

NonBond

R

Buckingham or exponential-6

E(R)= A e-CR - B/R6 which we write as

whereR0= Equilibrium bond distance;

r = R/R0 is the scaled distance

Do=well depth

z = dimensionless parameter related to force constant at R0

We define a dimensionless force constant as k= {(d2E/dr2)r=1}/D0

k = 72 for LJ12-6

z = 12 leads to –D0/r6 at long R, just as for LJ12-6

z = 13.772 leads to k = 72 just as for LJ12-6

A problem with exp-6 is that E- ∞ as R 0. To avoid this BioGraf, LinGraf, CeriusII calculate the inner maximum and reflect E about this point for smaller R so that E and E’ are continuous. When z>10 this point is well up the inner wall and not important

Converting between X-6 and LJ12-6

Usual convention: use the same R0 and D0 and set z = 12.

I recall that this leads to small systematic errors.

Second choice: require that the long range form of exp-6 be the same as for LJ12-6 (ie -2D0/r6) and require that the inner crossing point be the same. This leads to

This was published in the Dreiding paper but I do not know if it has ever been used

Popular vdW Nonbond Terms: Morse

NonBond

R

E(R) = D0 {exp[-b(R-R0)] - 2 exp[(-(b/2)(R-R0)]}

At R=R0, E(R0) = -D0

At R=∞, E(R0) = 0

D0 is the bond energy in Kcal/mol

R0 is the equilibrium bond distance

b is the Morse scaling parameter

I prefer to write this as

E(R) = D0 [X2 - 2X] where X = exp[(a/2)(1-r)]

At R=R0, r = 1 X = 1 E(R0) = -D0

At R=∞, X = 0 E(R0) = 0

k= {(d2E/dr2)r=1}/D0 = a2/2

Thus a = 12 k= 72 just as for LJ12-6

D0

R0

Few theorists believe that Morse makes sense for vdw parameters since it does not behave as 1/R6 as R∞

However, the vdw curve matches 1/R6 only for R>6A, and for our systems there will be other atoms in between. So Morse is ok.

vdW combination rules

Generally the vdw parameters are provided for HH, CC, NN etc (diagonal cases) and the off-diagonal terms are obtain using combination rules

D0IJ = Sqrt(D0II D0JJ)

R0IJ = Sqrt(R0II R0JJ)

zIJ = (zIJ + zIJ)/2

Sometimes an arithmetic combination rule is used for vdw radii

R0IJ = (R0II +R0JJ)/2 but this complicates vdw calculations

(the amber paper claims to do this but the code uses geometric combinations of the radii)

1-2 and 1-3 Nonbond exclusions

For valence force fields, it is assumed that the bond and angle quantities include already the Pauli Repulsion and electrostatic contributions so that we do not want to include them again in the nonbond list.

Thus normally we exclude from the vdw and coulomb sums

the contributions from bonds (1-2) and

next nearest neighbor (1-3) interactions

1-4 Nonbond exclusions

There is disagreement about 1-4 interactions. In fact the origin of the rotational barrier in ethane is probably all due to Pauli repulsion orthogonality effects.

Thus a proper description of the vdW interactions between 1-4 atoms should account for the barrier.

In fact it accounts only for 1/3 of the barrier.

I believe that this is because the CH bond pairs are centered at the bond midpoint not on the H atoms as assumed in the vdw.

Thus using atom centered vdw accounts for only part of the barrier necessitating an explicit dihedral term.

Thus I believe that the1-4 vdw and electrostatic terms should be included. However some FF, such as Amber, CHARRM reduce the 1-4 interactions by a factor of 2. I do not know of a justification for this.

Download Presentation

Connecting to Server..