Reaktivní silové pole ReaxFF
Tomáš Trnka
ttrnka@mail.muni.cz
Laboratoř výpočetní chemie (LCC)
Národní centrum pro výzkum biomolekul
Masarykova univerzita, Brno
C8855 2022-03-02
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 1 / 24
Silová pole, molekulová mechanika
▶ Klasická molekulová mechanika = “kuličky na pružinkách”
▶ Žádné explicitní elektrony, efektivní interakce atomů s okolím
▶ Pohyb atomů řízen klasickou mechanikou
▶ Potenciální energie systému je explicitní funkcí poloh všech atomů
▶ Potenciální energie = součet členů odpovídajících jednotlivým
interakcím / vnitřním stupňům volnosti
▶ vazby
▶ ohyb valenčních úhlů
▶ torzní úhly
▶ elektrostatické interakce
▶ repulze (elektrostatická repulze jader a elektronů, Pauliho repulze)
▶ disperzní interakce
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 2 / 24
Klasická silová pole
▶ Jednotlivé členy v MM potenciálu mají typicky jednoduchý empirický
tvar:
▶ vazby: harmonické pružiny Ebond,ij = k(rij − r0)
▶ elektrostatika: Coulombovská interakce dvou bodových nábojů
Eelst,ij =
qiqj
rij
▶ disperzní interakce: Edisp,ij = − C
r6
ij
▶ Nutno znát parametry (k, qi, C, …)
▶ Sada parametrů = “silové pole”
▶ Přenositelnost silových polí omezená na určitou skupinu systémů
▶ Nutnost zachytit různé chování atomů téhož prvku v závislosti na
chemickém okolí => více typů atomů s různými parametry pro jeden
prvek
▶ Mnoho kombinací (vazby: páry typů atomů, úhly: trojice, torze:
čtveřice) => stovky parametrů
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 3 / 24
Silová pole a chemické reakce
▶ Vazebné a úhlové členy jsou založeny na předdefinovaných vazbách
▶ Nelze popsat chemickou reakci (vazby nemůžou vznikat a zanikat)
▶ Typické řešení: hybridní metody (QM/MM)
▶ Reakční centrum popsáno kvantově-chemicky, zbytek MM silovým
polem
▶ Výpočetně náročné (nevhodné pro dlouhé MD simulace)
▶ Složité a náchylné k artefaktům
▶ Nevhodné pro velké systémy s mnoha reakcemi (max. stovky QM
atomů)
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 4 / 24
Reaktivní silová pole
▶ Empirické MM potenciály bez
předdefinovaných vazeb
▶ 1929: Morseův potenciál
▶ Popis vazebné vibrace diatomické molekuly
V(r) = D
(
1 − e−α(r−r0)
)2
▶ Parametry: vazebná energie D, rovnovážná
vzdálenost r0, tuhost pružiny kolem minima α
▶ Neomezená tvorba nových vazeb =>
nepoužitelný pro systémy obsahující více
atomů
0
0.5
1
1.5
2
2.5
3
3.5
4
0 1 2 3 4 5 6
U(r)/D
r/r₀
α = 1
α = 2
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 5 / 24
Reaktivní silová pole
▶ 1985 Abell, 1988 Tersoff: vazebné řády závislé na vzdálenostech
okolních atomů
▶ Síla vazebné interakce úměrná vazebnému řádu
▶ Součet vazebných řádů kolem atomu omezený jeho chemickou valencí
=> použitelné pro víceatomové systémy
▶ 1990 Brenner: REBO (Reactive Empirical Bond Order)
▶ Rozšíření Tersoffova potenciálu (pro Si) pro uhlovodíky
▶ První reaktivní FF pro organickou chemii
▶ 2000 Stuart: AIREBO (Adaptive Intermolecular REBO)
▶ REBO má špatně vdW disperzní interakce (Morseův potenciál není r−6
)
▶ Oprava: Kombinace REBO s Lennard-Jones 12-6
▶ LJ interakce kovalentně (1-2 a 1-3) sousedících atomů vyloučeny
(běžné v nereaktivních polích)
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 6 / 24
Silové pole ReaxFF
▶ 2001 Adri C. T. van Duin (Caltech, nyní Penn State):
reaktivní pole ReaxFF pro uhlovodíky
▶ Společné rysy s REBO, ale z většího odvozené od
univerzálního nereaktivního silového pole Delft
Molecular Mechanics (van Duin, 1994)
▶ 2008 Chenoweth et al.: finální tvar potenciálu ReaxFF
(FF pro oxidaci uhlovodíků – C/H/O)
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 7 / 24
Základní rysy ReaxFF
▶ Základ odvozený od konceptu centrálních sil (jen nevazebné interakce
všech párů) – dědictví DMM:
▶ Fluktuující bodové náboje počítané metodou EEM
▶ Morseův potenciál pro všechny páry atomů
▶ Dále doplněny lokální korekce pro správný popis větších molekul
▶ Vazby, úhly a dihedrály, vše úměrné vazebným řádům
▶ Žádné explicitní typy atomů - jen jedna sada parametrů pro každý
prvek
▶ Parametry fitovány na geometrie a křivky potenciální energie z QM
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 8 / 24
Členy potenciálu
U =Ebond + ELP + Eover + Eunder + Eval + Epen + Ecoa + EC2+
+Etriple + Etors + Econj + EH−bond + EvdWaals + ECoulomb
Ebond, Eval, Etors vazby/valenční úhly/dihedrální úhly
ELP energie volných elektronových párů
Eover, Eunder penalizace pro nad/pod-koordinaci každého atomu
Epen penalizace kumulovaných násobných vazeb (allen)
Ecoa třícentrová konjugace (karboxyly, nitroskupiny)
EC2 destabilizace biradikálů (C2, O2)
Econj čtyřcentrová konjugace a aromaticita
EvdWaals Morseův potenciál pro všechny páry
ECoulomb stíněná Coulombovská interakce EEM nábojů
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 9 / 24
Vazebné řády
BO′
ij = BOσ
ij + BOπ
ij + BOππ
ij = e
pbo1
(
rij
rσ
0
)pbo2
+ e
pbo3
(
rij
rπ
0
)pbo4
+ e
pbo5
(
rij
rππ
0
)pbo6
▶ “Surové” vazebné řády následně korigovány složitou procedurou tak,
aby atomy měly pokud možno očekávané valence a vazby chemicky
rozumné řády (0.2→0, 0.8→1, atd.)
▶ Energie vazby je přímo odvozena z korigovaných vazebných řádů (BO)
Ebond = −Dσ
e · BOσ
ij · epbe1(1−(BOσ
ij )
pbe2
) − Dπ
e · BOπ
ij − Dππ
e · BOππ
ij
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 10 / 24
Příklad: vazebný řád C-C
van Duin et al.; J. Phys. Chem. A, 2001, 105 (41), 9396–9409
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 11 / 24
Další členy potenciálu
▶ Úhlové a dihedrální členy pro každou dvojici/trojici sousedních vazeb
s BO > cutoff
▶ Rovnovážný valenční úhel závisí na sumě π vazebných řádů
(hybridizaci):
▶ Žádné π vazby: 109.47° (sp3), s rostoucím podílem π vazeb přes 120°
(sp2) až k 180° (sp)
▶ Torzní potenciál: kombinace kosinů s periodou 120° a 180°
▶ Druhý zmiňovaný příspěvek úměrný π charakteru centrální vazby
▶ Všechny nevazebné interakce hladce oříznuty (typicky na 10 Å)
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 12 / 24
Náboje v ReaxFF
▶ Interakce jen na krátkou vzdálenost (10 Å)
▶ Náboje atomů nemůžou být předdefinované (potřeba reagovat na
změnu konstituce a prostředí)
▶ Polarizovatelné silové pole – náboje počítány metodou EEM (Mortier
et al.: Electronegativity Equalization Method, též QEq - Charge
Equillibration)
EEEM =
∑
i
χiqi +
1
2
∑
i
ηiq2
i +
1
2
∑
i̸=j
qiqj
rij
▶ Dva parametry pro atom: elektronegativita χi a tvrdost ηi
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 13 / 24
Výpočet parciálních nábojů
▶ Minimalizace EEEM vázaná zachováním celkového náboje
▶ Soustava EEM lineárních rovnic
ηiqi +
∑
j̸=i
Jijqj = µ − χi Jij =
1
3
√
1
γ3
ij
+ r3
ij
▶ Interakce kovalentně vázaných atomů: nutnost stínění (Louwen-Vogt)
▶ Celkem tři parametry pro každý chemický prvek (χ, η, γ)
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 14 / 24
Struktura silového pole
▶ Žádné explicitní typy atomů - jen jedna sada parametrů pro každý
prvek
▶ Zhruba 30 parametrů pro každý chemický prvek, 20 pro každý pár
(vazbu), 5 pro úhel/torzi
▶ Vývoj silových polí složitý
▶ Přes 800 parametrů jen pro C/H/O/N
▶ Silné závislosti mezi jednotlivými členy
▶ Fitování na velké trénovací sady QM dat
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 15 / 24
Dostupné parametrizace
▶ Dvě základní větve (rodiny vzájemně částečně kompatibilních silových
polí)
▶ “Combustion branch”
▶ Reakce za vysokých teplot
▶ Reaktivita uhlovodíků, hoření, výbušniny, …
▶ “Water branch”
▶ Mírné podmínky, (nejen) vodné roztoky
▶ Degradace materiálů, polymery, …
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 16 / 24
Výhody ReaxFF
▶ Mnohem rychlejší než kvantově-chemické metody
▶ Zhruba 10-100x pomalejší než běžná nereaktivní silová pole
▶ 100k atomů, nanosekundy
▶ Jednoduchá příprava simulace
▶ Srovnatelné s QM, stačí prvky a souřadnice
▶ Netřeba přiřazovat typy atomů či definovat vazby
▶ Dostačující přesnost pro studium mnoha chemických procesů
▶ Široká škála dostupných silových polí pokrývajících většinu prvků
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 17 / 24
Omezení ReaxFF
▶ Špatný popis slabých nekovalentních interakcí
▶ Morseův potenciál místo správné disperze
▶ 10 Å cutoff i pro elektrostatiku
▶ Omezená spolehlivost EEM nábojů
▶ Častá nestabilita EEM modelu způsobující pády simulací
▶ Žádné explicitní elektrony => nelze popsat spinové stavy a redoxní
reakce
▶ Nejasná přenositelnost silových polí na nové systémy
▶ Náročná parametrizace a rozšiřování silových polí
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 18 / 24
Modifikace a rozšíření
▶ Nábojový model ACKS2 (T. Verstraelen, Gent)
▶ Složitější alternativa EEM řešící jeho hlavní problémy
▶ Neomezený přenos náboje (každý systém je vodivý)
▶ Neexistence chemicky rozumných iontových stavů (nesprávná disociace)
▶ Špatné škálování polarizovatelnost s velikostí systému
▶ Zatím málo rozšířený, většina FF používá EEM z historických důvodů
▶ ReaxFF-lg: disperzní korekce (Liu et al., W. Goddard; Caltech)
▶ Správný průběh disperzní interakce (r−6
)
▶ Zatím jen jedno silové pole s potřebnými parametry
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 19 / 24
Modifikace a rozšíření
▶ eReaxFF (Islam et al., PSU): pseudoklasické explicitní elektrony
▶ “Pseudoatomy” El/Eh pro elektrony a díry (bodový náboj ±1)
▶ Můžou se zcela volně pohybovat
▶ Interagují s “overcoordination/undercoordination” členy potenciálu
ReaxFF
▶ Studium organických radikálů aj.
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 20 / 24
Software pro ReaxFF
▶ Původní program “reaxff”
▶ Referenční implementace (van Duin)
▶ Pomalý, bez paralelizace, složité vstupy
▶ LAMMPS
▶ Svobodný SW
▶ Implementace ReaxFF převzata z PuReMD (Purdue, H. M. Aktulga)
▶ SCM: Amsterdam Modeling Suite
▶ Komerční
▶ Rychlejší než LAMMPS pro malé a střední systémy
▶ Podpora speciálních modifikací (eReaxFF, ACKS2)
▶ GUI, software pro vývoj silových polí (ParAMS)
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 21 / 24
ReaxFF a biomolekuly
▶ Parametrizace pro glycin (O. Rahaman) a peptidy/proteiny (S. Monti)
▶ Chybí popis správné dynamiky proteinů
▶ Typické aplikace: S. Monti – aminokyseliny/peptidy na zlatých
nanočásticích (drug delivery)
▶ LCC NCBR: vývoj metod pro studium enzymatických reakcí pomocí
ReaxFF
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 22 / 24
Vybrané aplikace ReaxFF
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 23 / 24
Cvičení: Spalování methanu
▶ $ cd /scratch/$LOGNAME
▶ $ module add adf
▶ $ export SCMLICENSE=/home/tootea/c8855.txt
▶ $ adfinput
▶ Překliknout oranžové políčko na “ReaxAMS”
▶ www.scm.com → Support → ReaxFF Tutorials → Reaction networks
▶ Následovat kroky 1-4
Tomáš Trnka (LCC NCBR MU) ReaxFF C8855 2022-03-02 24 / 24