Fyzika biopolymerů
Robert Vácha
Kamenice 5, A4 2.13
robert.vacha@mail.muni.cz
Základy simulací
"2
Force Field
google
"3
Force Field
- interakční potenciál … přibližný
- soubor rovnic a parametrů popisujících interakce mezi částicemi systému
- odvozen fitováním experimentálních a ab initio či jinak vypočtených dat
- označován za empirický potenciál (ale neni 100% protože obsahuje ab initio data)
- často dobře funguje jen na oblasti blízké fitovaným podmínkám (nejčastějí 300K, 1atm,
ve vodném roztoku, pH neutral, …)
- obvykle se takto označuje atomový model, lze chápat obecněji i na zhrubené modely
Properties	of	classical	mechanical	systems
5
• we	approximate	time	evolution	of	a	real	molecular	system	employing	 the		
classical	mechanics
• Newton’s	equations	of	motions	applied	to	all	atoms
• deterministic	time	evolution	(r0 and	v0 fully	determine	evolution)
• microscopic	time-reversibility
• macroscopic	time	irreversibility	as	an	emergent	behavior
• energy,	linear	momentum	 and	angular	momentum	conserved
• in	general,	chaotic	behavior	possible	
−
'(
')"
= $"
'*
)"
'+*
)" + = 0 , ."(+ = 0)
(()⃗) - force	field
)" + = 0 , ."(+ = 0)
)" + ,."(+) + = 0
+
Intermolecular	interactions
Intramolecular	interactions	classification:
• Electrostatic	interactions:	ion-ion,	 ion-dipole,	 ion-induced	dipole
• van	der	Waals	interactions	(all	but	the	above,	attractive		&	
repulsive,	between	uncharged	molecules):
• dipole-dipole	 (sometimes	classified	as	electrostatic)
• dipole-induced	 dipole	(Debye	force)
• between	two	instantaneously	induced	dipoles	(London	
dispersion	force)
Warning:	Various	classifications	exist
"4
Force Field all-atom
Vtot = Vbonded + Vnonbonded
Vbonded
Vnonbonded
( ) ( )
( )[ ] ∑∑∑
∑∑
+⎟
⎠
⎞
⎜
⎝
⎛
−+−+
+−+−=
elstat ij
ji
vdw ij
ij
ij
ij
dihed
n
anglesbonds
rtotal
R
qq
R
B
R
A
n
k
V
KrrKE
ε
γϕ
θθθ
612
cos1
2
0
2
0
Empirický potenciál
Vtot = Vbonded + Vnonbonded
"5
Force Field all-atom
Intramolecular	interactions	– bond	potential
Bond	potential	usually	approximated	as:
• Harmonic	potential
• Morse	potential
gromacs.org
BUT:	errors	due	to	differences	between	quantum	and	classical	oscillators!
IN	PRACITCE:	it	is	better	to	constraint	the	bonds	(it	also	allows	larger	Δt)
adapted	from	gromacs.org
8
Intramolecular	interactions	– bond	potential
Bond	potential	usually	approximated	as:
• Harmonic	potential
• Morse	potential
gromacs.org
BUT:	errors	due	to	differences	between	quantum	and	classical	oscillators!
IN	PRACITCE:	it	is	better	to	constraint	the	bonds	(it	also	allows	larger	Δt)
adapted	from	gromacs.org
8
- v praxi často drženo (constrained)
na pevné rovnovážné hodnotě pro
urychlení simulací
olecular	interactions	– valence	angle	potential
usually	approximated	as:
ngle	potential
o	as:
ed		potential	(in	GROMOS)
bending	potential
ey	potential	(in	CHARMM)
gromacs.org
9
Morseho potenciál
( )2)(
1) e
e
rra
eD −−
−=
potenciálové funkce:
arakterizuje šířku
ciálové jámy
loubka potenciálové jámy
ističtější reprezentace
y potenciální energie
harmonický potenciál
"6
Force Field all-atom
Intramolecular	interactions	– improper	dihedrals
Out-of-plane	bending,
used	to	maintain		planarity of	a	group
gromacs.org
11
Intramolecular	interactions	– torsion	(dihedral)	angle	potential
Different	(equivalent)	functional	forms:
• Periodic	(Fourier,	 proper)
• Ryckaert-Bellemans potential
gromacs.org
10
(2 3"456 = 7 89:
(1 + cos	(A3B − 3CB)
D
BEF
89:
- dihedral	force	constant
A – number	of	minimum
3CB - angular	offset
(usually	0ᵒ	for	odd	n,	180ᵒ for	even	n)
ramolecular	interactions	– torsion	(dihedral)	angle	potential
erent	(equivalent)	functional	forms:
Periodic	(Fourier,	 proper)
Ryckaert-Bellemans potential
gromacs.org
10
(2 3"456 = 7 89:
(1 + cos	(A3B − 3CB)
D
BEF
89:
- dihedral	force	constant
A – number	of	minimum
3CB - angular	offset
(usually	0ᵒ	for	odd	n,	180ᵒ for	even	n)
ntramolecular	interactions	– torsion	(dihedral)	angle	potential
Different	(equivalent)	functional	forms:
Periodic	(Fourier,	 proper)
Ryckaert-Bellemans potential
gromacs.org
(2 3"456 = 7 89:
(1 + cos	(A3B − 3CB)
D
BEF
89:
- dihedral	force	constant
A – number	of	minimum
3CB - angular	offset
(usually	0ᵒ	for	odd	n,	180ᵒ for	even	n)
Lennard-Jones	potential
Empirical	– fitted	to	experimental	data
another	form:
www.atomsinmotion.com/book
"7
Force Field all-atom
Electrostatic	interactions	– charge-charge
The	longest	range!
("4 =
G"G4
4IJF)"4
r
15
Lennard-Jones	potential
Empirical	– fitted	to	experimental	data
another	form:
www.atomsinmotion.com/book
"8
Potential energy surface
Born-Oppenheimerova aproximace: elektrony se okamžitě přizpůsobují pohybu jader
(jádra jsou mnohem těžší než elektrony),pohyb jader tedy lze řešit klasicky a pro každou
konfiguraci jader vyřešit Schrodingerovu rovnici = vlnovou funkci rozložení elektronů
Born-Oppenheimerova aproximace
Separace pohybu elektronů a pohybu jader
Jako nezávislé proměnné hamiltoniánu jsou pouze elektrony, pozice jader
jsou považovány za fixované (jádra jsou mnohem těžší než elektrony, na
škále elektronových vibrací se poloha jader mění velmi pomalu)
Koncept hyperplochy potenciální energie (potential energy surface
PES) = hyperplocha potenciální energie přes všechna možná
uspořádání atomů
Dimenze PES: 3N-6
Mnoharozměrný (počet stupňů
volnosti) povrch potenciální energie
po kterém se částice pohybují ,,,
může být popsán pomocí kvantové
teorie, forcefiledu, modelu…
Pro řešení pohybu v kvantové teorii separujeme pohyb elektronů a jader
"9
Molekulová dynamika (MD)
- řeší klasické Newtonovy rovnice
- výsledkem je trajktorie = set po sobě jdoucích konfigurací systému s definovaným
časovým odstupem
- vlastnosti systému se počítají analýzou trajektorie
- musí být zadány počáteční podmínky - polohy a rychlosti všech částic
- používá forcefield, který musí mít definované síly (první derivace potenciálu)
- deterministická
- časově reversibilní
- zachování momentů hybnosti
- numerické řešení rovnic pomocí integrátoru
Key	points	from	the	previous	lecture
lassical	Newton’s	
ns	of	motion	are	solved	
ms	in	the	considered	
ry	(positions	of	atoms	
is	the	result	of	MD	
on.	System	properties	
ed	from	the	trajectory.
nditions	must	be	set	
ositions	and	velocities	
oms).
on	potential	must	be	
d	(force	field).
!" = $"%"
−
'(
')"
= $"
'*
)"
'+*
)" + = 0 , ."(+ = 0)
(()⃗) - force	field
Properties	of	classical	mechanical	systems
5
• we	approximate	time	evolution	of	a	real	molecular	system	employing	 the		
classical	mechanics
• Newton’s	equations	of	motions	applied	to	all	atoms
• deterministic	time	evolution	(r0 and	v0 fully	determine	evolution)
• microscopic	time-reversibility
• macroscopic	time	irreversibility	as	an	emergent	behavior
• energy,	linear	momentum	 and	angular	momentum	conserved
• in	general,	chaotic	behavior	possible	
−
'(
')"
= $"
'*
)"
'+*
)" + = 0 , ."(+ = 0)
(()⃗) - force	field
)" + = 0 , ."(+ = 0)
)" + ,."(+) + = 0
+
Properties	of	classical	mechanical	systems
5
• we	approximate	time	evolution	of	a	real	molecular	system	employing	 the		
classical	mechanics
• Newton’s	equations	of	motions	applied	to	all	atoms
• deterministic	time	evolution	(r0 and	v0 fully	determine	evolution)
• microscopic	time-reversibility
• macroscopic	time	irreversibility	as	an	emergent	behavior
• energy,	linear	momentum	 and	angular	momentum	conserved
• in	general,	chaotic	behavior	possible	
−
'(
')"
= $"
'*
)"
'+*
)" + = 0 , ."(+ = 0)
(()⃗) - force	field
)" + = 0 , ."(+ = 0)
)" + ,."(+) + = 0
+
"10
Molekulová dynamika (MD)
- Historie
1957 - Alder and Wainwright - basics of MD
1964 - Rahman - MD with Lennard-Jones potential, NVE, liquid argon
1967 - Verlet - Verlet integration algorithm and Verlet neighbor list
1974 - Stillinger and Rahman - MD of water
1981 - Andersen and Parrinello-Rahman - NPT
1986 - Nose and Hoover - NVT with Nose-Hoover thermostat
1985 - Car and Parrinello - ab initio MD
"11
Integrátory = numerické řešení pohybových rovnic
Taylorův rozvoj
Euler	algorithm
8
NOT	USED	in	MD:
• time	irreversible	(but	Newton’s	equations	are	reversible)
• no	phase-space	preserving	(Liouville’s theorem	violated)
(forward	difference	approximation)
není časově reverzibilní a nezachovává fázový objem
Time	irreversibility	of	integrators
9
in	forward	
propagation	forces	
calculated	at	+
+ + ∆++
) + = 0 , .(+ = 0)
) + ,.(+)
∆+
in	backward	
propagation	forces	
calculated	at	+ + ∆+
as	a	trick	to	symmetrize	
integrator,	forces	can	be	
calculated	in	the	middle	
+ → −+;. + → −.(+)
Luivillův teorém
4
Liouville’s	theorem
4
'E(=,4)
'+
=
2E
2+
+ F
2E
24"
4"̇ +
2E
2="
="̇
8
"GH
= 0
4
=
+ = 0 +
∆="∆4" = JKLM+note:
The	distribution	 function	 is	constant	along	any	trajectory	in	phase	space.
"12
Integrátory
Verlet algorithm
10
Derivation:
Positions	from	adding	and	velocities	from	subtracting	them:
coordinates
+ + ∆++ − ∆+
new	coordinatescoordinates
&	forces
+
Verlet alg.	is	time	reversible	(symmetric	when	+ → −+)
Verletův algoritmus
Verlet algorithm
10
Derivation:
Positions	from	adding	and	velocities	from	subtracting	them:
coordinates
+ + ∆++ − ∆+
new	coordinatescoordinates
&	forces
+
Verlet alg.	is	time	reversible	(symmetric	when	+ → −+)
řešením reverzibility je symetrizace
sečtením rovnic….
Velocity-Verlet integrator
14
(available	in	Gromacs)
• positions,	 velocities,	and	forces	available	at	the	same	time
• time-reversible	and	phase	space	volume	preserving
• 2nd order	but	more	accurate	than	Verlet,	very	stable
• better	numerical	accuracy	because	velocities	are	stored
• in	practice,	an	optimal	integrator	for	MD
• (always	check	your	MD	software	recommendations!)	 "13
Integrátory
Potřeba termostatu, barostatu, atd. !!! specifické přednášky
Rychlostní-Verletův algoritmus (Velocity-Verlet)
- je časově reverzibilní
- zachovává fázový objem => ideální řešení
- velmi stabilní ale zkontrolujte doporučení
používaného programu
Leap-frog integrator
13
(default	in	Gromacs)
• velocities	updated	at	half	time	steps	and	‘leap’	ahead	the	positions
• 3rd order	algorithm,	time	reversible,	phase-space	volume	not	preserved
• algebraically	equivalent	to	Verlet
• current	velocities	can	be	approximated	as:
Leap-frog (Skákající žáby) algoritmus
13
• velocities	updated	at	half	time	steps	and	‘leap’	ahead	the	positions
• 3rd order	algorithm,	time	reversible,	phase-space	volume	not	preserved
• algebraically	equivalent	to	Verlet
• current	velocities	can	be	approximated	as:
algebraicky stejné řešení jako Verlet
- časově reverzibilní
- ale nezachovává fázový objem
"14
Příklady Molekulové dynamiky
"15
Příklady Molekulové dynamiky
"16
Příklady Molekulové dynamiky
"17
Příklady Molekulové dynamiky
"18
Monte Carlo (MC)
- obecná metoda na vzorkovaání pomocí náhodných čísel
- molekulární x nemolekulární (výpočet integrálů, diferenciálních rovnic atd.) používá
se ve finančnictví, počítačové grafice nebo umělé inteligenci
Historie
• 1930s & 1940s – Fermi, Ulam, von Neumann: Los Alamos, the Manhattan
project = procházení neutronů skrz materiál
• 1950s – Teller: vývoj vodíkové bomby
• 1953 – Metropolis (Metropolis-Hastings algorithm) “Equation of State
Calculations by Fast Computing Machines”
• 1966 – Whittington et al. “Effect of density on configurational properties of longchain
molecules using a Monte Carlo method”
• 1987 – Li & Scheraga: Monte Carlo - minimalizace pro hledání mnoha minim v
proteinovém sbalování (foldingu)
• 1992 – Leach: Dockování ligandu do proteinů s flexibilitou postraních řetězců
• 2011 – Chodera: Pohyb nerovnovážného kandidáta
• 2011 – Martinez Veracoechea & Frenkel : Hybridizační pohyb vytvoření vazby
mezi ligandy a receptory
"19
Příklad
Vypočtěte pi pomocí nahodného vzorkování bodu ve čtverci a jemu vepsanému kruhu
"20
Řešení
Monte	Carlo	– calculation	of	circle	surface	area
5
2r
!"#$%&' = (2+)!./&.0'
=	?
Computation:
1. Generate	3 random	points	in	the	square
2. Calculate	number	of	points	in	(4/5)	the	circle
Probability	to	hit	the	circle	6/5 =
789:;8<
7=>?@:<
=	
59A
B
Hence:		!./&.0' = !"#
59A
B
= 4+- 59A
B
(valid	for	3 → ∞)
Similar	algorithms	can	be	used	for	calculation	of,	for	example:
• integrals
• surface	areas	of	object	with	complex	shapes	(e.g.,	proteins)
• volumes	of	protein	channels
Pravděpodobnost, že náhodný pod je v kruhu
P =
Akruh
Actverec
P =
R2
(2R)2
= /4
Obdobně lze spočítat objem velmi komplikovaných objektů
"21
Kanonický soubor
- vzorkování konfigurací pomocí nahodných změn
- pravděpodobnost konfigurací je dán Boltzmannovým vztahem
P(E) e
E
kT
- střední hodnota veličiny
Monte	Carlo	– canonical	ensemble	(direct	sampling)
6
Computation	of	any	macroscopic	quantity:
1. Generate	a	series	of	microstates	(with	different	ri,	vi)	with	
constant	N,	V,	T	(these	states	represent	in	the	same	macrostate)
2. Probability	of	finding	 a	configuration	with	energy	E	is	
propotional	 to	the	Boltzman	factor	(from	canonical	ensemble	
properties):
6 F ~HIJ −
F
LM
3. Mean	value	of	any	macroscopic	quantity:
! =
∑ !/HIJ −
F/
LM
∑ HIJ −
F/
LM
But	where	are	the	random	numbers?
The	different	microstates	can	be	generated	
by	random	changes	of	positions	and	
velocities	of	particles
Monte	Carlo	– canonical	ensemble	(direct	sampling)
7
Problem:
For	randomly	generated	microstates,	many	configurations	lead	to	
atomic	overlaps	thus	giving	high	energy	and	low	Boltzman	factor	
and	hence	not	contributing	 much	to	sums	in	the	average:
! =
∑ !/HIJ −
F/
LM
∑ HIJ −
F/
LM
Becuase,	in	other	words,	in	‘real’	canonical	enseble:
6 F ~HIJ −
F
LM
Solution:
Find	a	way	to	generate	‘important’	confgurations	 and	avoid	‘not	
important’	ones
P
E
Monte	Carlo	– canonical	ensemble	(direct	sampling)
7
Problem:
For	randomly	generated	microstates,	many	configurations	lead	to	
atomic	overlaps	thus	giving	high	energy	and	low	Boltzman	factor	
and	hence	not	contributing	 much	to	sums	in	the	average:
! =
∑ !/HIJ −
F/
LM
∑ HIJ −
F/
LM
Becuase,	in	other	words,	in	‘real’	canonical	enseble:
6 F ~HIJ −
F
LM
Solution:
Find	a	way	to	generate	‘important’	confgurations	 and	avoid	‘not	
important’	ones
P
E
- problém se vzorkováním - většinu času se vzorkují stavy, které přispívají
minimálně, důležité stavy jsou málo vzorkovány
"22
8
E ! =
∑ !/HIJ −
F/
LM
∑HIJ −
F/
LM
P
E
Importace	(weighted)	sampling
Sample	more	often	in	the	areas	with	
low	energy	contributing	 more	to	the	
sum!
Metropolis-Hastings algorithm
- vzorkuje se oblast (kofigurace), které jsou v kanonickém souboru důležité
- nové konfigurace jsou příjmuty s pravděpodobností
There are other formulas that satisfy this relation, but the standard Metropolis
condition is:
p(old ! new)
= P(new)
P(old)
if P(new) < P(old)
= 1 if P(new) P(old)
(10.5)
We will demonstrate this method on a simple system. Assume that we
have a system in a NV T ensemble with two particles A, B moving on the
x-axis with simple a displacement move. The displacement is ±1 with the
same probability in both directions. The configuration probability is then
P ⇠ exp( E/kT) and the probabilities of the displacement are:
p(old ! new) = e
Enew Eold
kT if Enew > Eold
= 1 if Enew  Eold
(10.6)
Příklad algoritmu:
1. náhodně vyber částici
2. vyber náhodný směr
3. pohni částicí
4. příjmi novou polohu na základe pravděpodobnosti výše
"23
Jak navrhnout nový krok
- Podmínka detailní rovnováhy (detailed balance condition)
94
Let’s start with basic Metropolis MC for canonical ensemble. The sampling
starts from an initial configuration with energy Eold. A new configuration
with energy Enew is generated by a move (for instance a small displacement)
and this new configuration is accepted or rejected based on the
transition probability. There are several options for calculating the transition
probability depending on the particular move, but in general the transition
probability has to keep the system in equilibrium once it is reached. In other
words, once the system is in equilibrium the number of visits to particular
configurations is proportional to the configuration probability. We usually
use the much stronger condition of detailed balance, where the averaged
number of moves from the old state to the new one is equal to the average
number of reverse moves. This can be formulated as:
P(old)a(old ! new) = P(new)a(new ! old) (10.2)
where a is the transition probability and P the configuration probability. The
ratio of transition probabilities is then:
a(old ! new)
a(new ! old)
=
P(new)
P(old)
(10.3)
Each transition probability a consists of the probability of generating the
given move ↵ and the probability of accepting such a move p. In most of
the common methods ↵(old ! new) = ↵(new ! old) and hence we can
write the acceptance ratio:
p(old ! new)
p(new ! old)
=
P(new)
P(old)
(10.4)
There are other formulas that satisfy this relation, but the standard Metropolis
condition is:
p(old ! new)
= P(new)
P(old)
if P(new) < P(old)
= 1 if P(new) P(old)
(10.5)
We will demonstrate this method on a simple system. Assume that we
have a system in a NV T ensemble with two particles A, B moving on the
x-axis with simple a displacement move. The displacement is ±1 with the
same probability in both directions. The configuration probability is then
P ⇠ exp( E/kT) and the probabilities of the displacement are:
p(old ! new) = e
Enew Eold
kT if Enew > Eold
= 1 if Enew  Eold
(10.6)
P je pravděpodobnost stavu, a je pravěpodobnost přechodu skládající se z
pravděpodobnosti výběru daného pohybu krát pravděpodobnost přijetí daného
pohybu p
v rovnováze se
zastoupení stavů nemění
94
Let’s start with basic Metropolis MC for canonical ensemble. The sampling
starts from an initial configuration with energy Eold. A new configuration
with energy Enew is generated by a move (for instance a small displacement)
and this new configuration is accepted or rejected based on the
transition probability. There are several options for calculating the transition
probability depending on the particular move, but in general the transition
probability has to keep the system in equilibrium once it is reached. In other
words, once the system is in equilibrium the number of visits to particular
configurations is proportional to the configuration probability. We usually
use the much stronger condition of detailed balance, where the averaged
number of moves from the old state to the new one is equal to the average
number of reverse moves. This can be formulated as:
P(old)a(old ! new) = P(new)a(new ! old) (10.2)
where a is the transition probability and P the configuration probability. The
ratio of transition probabilities is then:
a(old ! new)
a(new ! old)
=
P(new)
P(old)
(10.3)
Each transition probability a consists of the probability of generating the
given move ↵ and the probability of accepting such a move p. In most of
the common methods ↵(old ! new) = ↵(new ! old) and hence we can
write the acceptance ratio:
p(old ! new)
p(new ! old)
=
P(new)
P(old)
(10.4)
There are other formulas that satisfy this relation, but the standard Metropolis
condition is:
p(old ! new)
= P(new)
P(old)
if P(new) < P(old)
= 1 if P(new) P(old)
(10.5)
We will demonstrate this method on a simple system. Assume that we
have a system in a NV T ensemble with two particles A, B moving on the
x-axis with simple a displacement move. The displacement is ±1 with the
same probability in both directions. The configuration probability is then
P ⇠ exp( E/kT) and the probabilities of the displacement are:
p(old ! new) = e
Enew Eold
kT if Enew > Eold
= 1 if Enew  Eold
(10.6)
- pokud je pravděpodobnost výběru daného pohybu vpřed a zpěz stejná existuje
jednoduché pravidlo splnující podmínku rovnováhy (mohou být i další)
např. v NPT (isobaricko-isochorickém) souboru
96
and the MC method samples this integral. In the same way we will use an
MC scheme to sample the integral of the partition function of the isobaric
ensemble, which is:
Z(N, p, T) =
Z
Q(N, V, T)e
pV
kT dV
=
p
⇤3N N!kT
Z
V N
e
pV
kT e
E
kT drN
dV (10.11)
In other words the density of states (the probability that a system of N
particles will have volume V at pressure p) is:
P(N, V ) = V N
e
pV
kT e
E
kT (10.12)
= exp(
E + pV NkT ln(V )
kT
) (10.13)
and the probability criterion for accepting or rejecting the volume change
is:
prob(old ! new) = min
⇢
1, exp

Enew Eold
kT
+
p(Vnew Vold) NkT ln(Vnew/Vold)
kT
(10.14)
We use prob for probability here to avoid the confusion with pressure p.
Usually it is more efficient to make random changes in the logarithm of
volume (ln(V )) than in the volume itself. The acceptance criterion then
has to be modified to:
prob(old ! new) = min
⇢
1, exp

Enew Eold
kT
+
p(Vnew Vold) (N + 1)kT ln(Vnew/Vold)
kT
(10.15)
since
R
V N
dV =
R
V N+1
d(ln(V )).
If we want to study adsorption in computer simulations, we can either
perform a simulation in the canonical ensemble with a system large enough
to represent infinite bulk, or we can employ a Grand Canonical ensemble
simulation. Simulation with the Grand Canonical ensemble only needs to
simulate the surface, therefore the system is much smaller and faster to calculate.
The bulk reservoir with constant temperature and pressure (chemical
potential) is replaced by a fluctuating number of particles in the system.
p(old new) = min 1, exp
Enew Eold
kT
+
p(Vnew Vold) NkT ln(Vnew/Vold)
kT
"24
Kinetické Monte Carlo
Kinetic	Monte	Carlo
14
Kinetic	Monte	Carlo
Assumption:
the	rate	of	each	process	depends	on	its	activation	
energy	(barrier):
+/ ∝ exp	 −
F`/
LM
Algorithm:
1. Being	in	the	current	state,	chose	randomly	one	
process	possible	to	occur	from	this	state
2. Realize	the	process	with	probability:
6/ =
+/
∑ +a
Initial	state
r1
r2
r3
r5r4 r6
r7
Advantages:
non-equilibrium	 and	non-stationary	processes	
can	be	simulated!
- musíme znát rychlostní konstanty např.
Kinetic	Monte	Carlo
14
Kinetic	Monte	Carlo
Assumption:
the	rate	of	each	process	depends	on	its	activation	
energy	(barrier):
+/ ∝ exp	 −
F`/
LM
Algorithm:
1. Being	in	the	current	state,	chose	randomly	one	
process	possible	to	occur	from	this	state
2. Realize	the	process	with	probability:
6/ =
+/
∑ +a
Initial	state
r1
r2
r3
r5r4 r6
r7
Advantages:
non-equilibrium	 and	non-stationary	processes	
can	be	simulated!
e
EB
kT
Algoritmus:
1. náhodně se vybere proces
2. proces se zrealizuje s pravděpodobností
3. aktualizuj čas
- existuje algoritmus bez zamítání procesů
- lze simulovat nerovnovážné procesy
Dynamické Monte Carlo
- musíme znát difuzní koeficienty (translační, rotační,..)
- parametry pohybů (maximalní velikost pohybu, pravdepodobnost pohybu) se
nastaví podle difuznich koeficientů
- konverguje k difusivnímu (brownovskému) pohybu
"25
Příklady Monte CarlaExamples	of	Monte	Carlo	methods	I
10
Ising model (originally	 for	spins)
Ising model	for	mixed	lipid	membranes
Lis	et	al.	JCTC	2012
"26
Příklady Monte Carla
Examples	of	Monte	Carlo	methods	II
11
Polymer	on	a	lattice
Protein	folding
Polymer	folding
"27
Příklady Monte Carla
"28
Příklady Monte Carla
"29
MD x MC
MD
výhody:
- dynamika vývoje systému
automatickou součástí
- informace o rychlostech
momentu atd
- rovnovážná i nerovnovážná
- jednoduchý algoritmus
- spousta optimalizovaného
softwaru
nevýhody:
- pouze malé změny v
jednotlových krocich = pomalost
- nutnost sil = pouze spojitě
diferencovatelné potenciály
- potřeba dodatečných algoritmů a
nastavení parametrů pro jiné než
mikrokaonické soubory
MC
výhody:
- obecnost a variabilita (systém,
potenciály, termodynamický soubor,
pohyby a změny systému)
- nevyžaduje síly
- možnost přeskakovat bariéry, snadné
úniky minim
- algoritmus více flexibilní
- snadno různé termodynamické
soubory
nevýhody:
- systémy v rovnováze a bez dynamiky
(nerovnováha a dynamika vyžadují
speciální přístup)
- programy nejsou optimalizované,
často vyžadují vlastní úpravu
"30
Ergodická hypotéza
= průměr získaný pro malý počet molekul přes dlouhý čas je ekvivalentní
průměrování přes velký počet molekul a krátký čas (v limitě: časový průměr pro
jednu molekulu je ekvivalentní průměru pro velký počet molekul)
- výsledky simulací a experimentů představují průměrné hodnoty (průměrování
přes počet molekul, čas)
Ergodická hypotéza
• výsledky simulací a experimentů představují průměrné
hodnoty (průměrování přes počet molekul, čas)
• např. IR spektrum: 1018 molekul, časový úsek 10-14 s
simulace: 103 molekul, časový úsek 10-9 s
• ergodická hypotéza: průměr získaný pro malý počet
molekul přes dlouhý čas je ekvivalentní průměrování
přes velký počet molekul a krátký čas (v limitě: časový
průměr pro jednu molekulu je ekvivalentní průměru pro
velký počet molekul)
∑∫ =
∞→∞→
==
M
i
i
M
X
M
dttXX
10
1
lim)(
1
lim
τ
τ τ