Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.elch.chem.msu.ru/mole/lecture/probst.pdf Äàòà èçìåíåíèÿ: Tue Aug 28 21:47:04 2012 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 23:01:01 2012 Êîäèðîâêà:

MD Simulations: Introduction and Ethylene Glycol as Example
·

Michael Probst1, Renat Nazmutdinov2, Oksana Ismailova1, Alexander Kaiser1
· · · Inorganic Chemistry Kazan State Technological University Kazan,Tatarstan, Russian Federation

· Institute for Ion Physics and Applied Physics · University of Innsbruck, Austria

2012-08-28 Dubna


Why does one need MD ­ Simulations ?
·

The interactions within a chemical system (Atoms, molecules, ions ...) result basically from its electronic structure which governs the interaction between the particles. Suppose, the total energy E(x) is known as a function of the coordinates x of all particles: How can one calculate thermodynmamic properties from E ? Only by numerical simulation.

·

·

·


Why numerical ?
· Reason 1: Integration over phase space:
· Integral over a high-dimensional weighted function

· Reason 2: Calculation of trajectories.
· analytically solvable only for E(x)=k x2 :

X(t)=a sin(wt)+b cos(wt)


MD principle 2

Many variations:
Accurate Slow on the computer Large Systems Simple Fast Small Systems

Accurate
Exact QM

Simple
(Semi) Classical Mechanics Approximate QM Interaction potentials electron densities Continuum Point Charges

Electrons Atoms


CC principle 4

CC & Molecular Simulations


MD principle 3

Other __________________________________ are very much related to MD


MD priciples

Advantage: One can calculate large systems (thousands and millions of atoms). Disadvantage: The potential energy function must be provided. It`s quality is crucial for the result. Bridging the time gap (like the spacial gap before):
Quantumchemical Hamiltonian Analytical Potential Function Characteristic Diff.equ.

100 particles ps simulation

105 particles ns simulation

mesoscale etc.


EG

Example:

Properties of Ethylene glycol in the condensed phase

· Flexible and symmetric H-O-C-C-O-H chain structure · Molecular physical properties of ethylene glycol depend on chain angles and aggregation


A)

Interaction energy

bond stretch

torsional

intermolecular interactions valence angle bend

intramolecular nonbonded


A bit more detail
Intramolecular: ·Bonded ·Non-bonded

E

bonded

=E

bond - stretch

+E

angle -bend

+E

rotate - along -bond


E

bond - stretch

=

1, 2 pairs



K b (b - b0 )

2

E

bond - bend

=

angles



K ( - 0 )

2

E

rotate - bond

=

1, 4 pairs



K (1 - cos(n ))


Interactions between non-Bonded Atoms
van der Waals · Electrostatic
·

E

non- bonded

=E

van - der -Waals

+E

electrostatic

Van der Waals: A Lennard-Jones (LJ) form is a compromise. Powers can vary.

Electrostatic / Coulomb: Partial charges on atoms.

E

Lennard - Jones

=

nonbonded pairs



Aik C ik 12 - 6 r rik ik

E

electrosta tic

=

nonbonded pairs



qi qk Drik


For EG:


Simulation = Trajectory calculation · We need the movements (=trajectories x(t))

B)

of the atoms and have to solve (`integrate') their equations of motion. · That means to calculate x(t) from

·

Since E(x) is complicated, this can only be done numerically which is, however, quite easy.


The equations of motion
Many possibilities: · Verlet: r(t+t) = 2r(t) - r(t-t) + t2a(t) · Leapfrog: r(t+t) = r(t) + t v(t+t/2)v(t+t/2) = v(t-t/2) + t a(t) · Gear ...
· ·

Velocity verlet algorithm: r(t+t) = r(t) + t v(t) + t2a(t) v(t+t) = v(t) + t [a(t) + a(t+t)]/2


Equilibration :

NVE ensemble at 298 K. 512 EG molecules.


Typical protocol:
· 512 EG molecule in cubic periodic box.

· 1) Equilibration in NVE ensemble at 298K for 5 ps.
· Time step of 0.5 fs.

· 2) Equilibration in NVT ensemble at 298K for 5 ps with a
weak-coupling Nose-Hoover thermostat(=modified EOMs). · Simulation time 500 ps.

· 3) Production run in NVT at 298K for 4 ns.
· Time step of 1 fs. · Trajectory recording frequency of 2.


C)

Analysis
Analysing the trajectories

Thermodynamic averages

`normal' averages ­ just statistics:

A=

lim
t

1 A(t )dt = t

lim
m i =1

1 m

m

Ai

Average over timesteps (and within a timestep)


Thermodynamic averages
Analysing the trajectories Example:

Average value of dipole moment T (K) Dipole moment (D) Experiment
Other example: Diffusion constant

200 250 298 350 400 2.40 2.37 2.48 2.59 2.80 2.25


The Distribution functions:

radial distribution function

1 g (r) = N


t =1 j i

N

TS

N

pair

( r - rij )


The Distribution functions:

radial distribution function


CC def 2

Work mode ­ computer related
1)
· Perform simulation, create trajectory file · Calculate properties timestep by timestep

2)
· Perform simulation, create trajectory file · Copy it all (as some kind of object, possibly in chunks) to memory · Calculate properties by filtering operations.


CC def 3

Angular distributions are similar:

O-C-C-O

Distribution of the O-C-C-O dihedral angle.


CC def 2

Dynamical properties:
·

Properties as a function of time.
H6-O3-C2-C1


NPT ENSEMBLE AT 200K O-C-C-O

a) 1.4% H-O-C-C

T

b) 0.8%

c) 4.3%

d) 4.5%


NPT ENSEMBLE AT 298K O-C-C-O

a) 20.3% H-O-C-C

T

b) 16.8%

c) 11.7%

d) 10.1%


NPT ENSEMBLE AT 400K O-C-C-O

a) 25.0% H-O-C-C

T

b) 21.7%

c) 15.9%

d) 17.7%


CC def 3

·

General idea: Standard procedure:
· A continuous distribution is discretized into (G and T)

T

· This frees up one dimension for analysing / visualizing it (`Histogram')
· Can be combined with another tool: nD histograms:


CC def 3

· DISTRIBUTION OF THE (static) CORRELATION BETWEEN THE 2 H-O-C-C

DIHEDRAL ANGLES:

NPT ENSEMBLE AT 400K

T
· For O-C-C-O in G, the probability of finding one conformation is slightly larger. · Both in trans are unlikely. H-O-C-C in trans


CC def 3

· DISTRIBUTION OF THE (static) CORRELATION BETWEEN THE 2 H-O-C-C

DIHEDRAL ANGLES:

NPT ENSEMBLE AT 400K

T
· For O-C-C-O in G, the probability of finding one conformation is slightly larger. · Both in trans are unlikely. H-O-C-C in trans


CC def 3

·

D) Visualisation of MD results
· Serious science: Extracting information of high-dimensional data · One good program: VMD


CC def 3

In our EG ­ work: · Extremely difficult problem: Undrstanding the H-bond network over space and time. · Our own program: VISH (W. Benger et al.)


D)
·

The interesting issues in electrochemistry
Solvated molecules that can have different oxidation states S2O82- + e S2O83SO42- + SO4-

· ·

Interesting are the reorganisation energies:

s = En

oneq

( S 2O ) - Eeq ( S 2O )

2- 8

2- 8

s = Enoneq ( S 2O83- ) - Eeq ( S 2O83- )


Oxidation / redaction reactions
·

From the RDFs of S2O8 in EG / H2O ...

Coordination number Distance of first minimum of RDF (nm)

water S2O82- S2O8325 0.58 26 0.56

EG S2O8213 0.673 S2O8312 0.638

·

... one can calculate the solvent reorganisation energies:
x(EG)
s (red) s (ox) s (avg)

0 1.0

22.4 38.3

kcal mol 54.5 32.2

-1

33.2 34.9

·

The increase of the average s for EG means that a dynamic solvent effect modulates the ET.


Thank you ...