Молекулска динамика (МД) је облик рачунарске симулације, где атоми и молекули могу да интерагују у одређеном временском интервалу, покоравајући се познатим законима физике. Молекулски системи се генерално састоје од великог броја честица^[1], стога је немогуће аналитички пронаћи особине тако сложених система. Молекулска динамика превазилази овај проблем користећи нумеричке методе. Чинећи везу између лабораторијских експеримената и теорије, можемо је посматрати као виртуелни експеримент.

Иако је Болцман (енгл. Ludwig Boltzmann) још у 19. веку објаснио да се материја састоји од интерагујућих честица које се крећу, већина људи и даље молекуле сматра ригидним. Ричард Фајнман (енгл. Richard Feynman) је 1963. рекао: „Све што жива бића раде, можемо схватити као играње и скакутање атома“.^[2] Једно од разјашњења које је молекулска динамика начинила, је развитак свести да су молекули, као што су протеини и ДНК, покретни системи.^[3] Молекулска динамика објашњава везу између молекулске структуре, њихове функције и кретања молекула.

МД је мултидисциплинарног карактера. Базира се на примени математичких алгоритама и информационе теорије (дисциплина примењене математике). Настала је у оквиру теоријске физике 1950. године, док се данас углавном примењује у науци материјала и у раду са биомолекулима.

Пре него што је омогућена рачунарска симулација молекулске динамике, овај задатак је решаван применом физичких модела, као што је макроскопска сфера. Идеја је била да се молекули сложе тако да приближно описују особине течности. 1962. Џон Дезмонд Бернал (енгл. John Desmond Bernal) је рекао: ”... узео сам неколико гумених лопти, сабио их тако да је међусобно растојање износило од 2,75 инча до 4 инча. Покушао сам то да урадим што смиреније, радећи у својој канцеларији, али на сваких пет минута био сам прекидан тако да нисам успео да запамтим шта сам пре тога радио”.^[4] На срећу, данас компјутери чувају информације о везама током симулације.

Молекулска динамика је специјална дисциплина молекулског моделирања и рачунарске симулације, базирана на статистичкој механици. Главни разлог коришћења МД методе је тај што је просечан статистички ансамбл једнак просечном времену система, познато као ергодска хипотеза. МД је такође означена као "статистичка механика бројева" и "Лапласова визија Њутнове механике" предвиђања будућности помоћу анимације, омогућава поглед на кретање атома унутар молекула. Са друге стране, дужа МД симулација је математички лоше условљена. Начињене грешке у нумеричкој интеграцији могу бити минимизиране одговарајућим одабиром алгоритама и параметара, али не могу бити потпуно уклоњене. Постојеће потенцијалне функције су у многим случајевима недовољно тачне за приказ динамике молекулских система. Како год, МД омогућава увид у прецизно време и просторни распоред унутар фазног простора, простора унутар кога су могућа сва стања система.

Постоји значајна разлика између области и метода које користе хемичари и физичари, и то се огледа у различитој терминологији коју користе. У хемији и биофизици, интеракција између честица се преноси "пољем сила" (класична МД), квантним хемијским моделом или комбинацијом ова две могућности. Ови термини се не користе у физици, где се интеракција описује под називом теорија или апроксимација коју називамо потенцијална енергија или само потенцијал.

Са почетком у теоријској физици, метод МД је постао популаран у науци материјала, а од 1970-их такође и у биохемији и биофизици. У хемији има кључну улогу у одређивању структуре протеина, и при употреби метода као што су X-зраци у кристалографији и НМР. Такође се ограничено употребљава у методи предвиђања структуре протеина. Користи се у испитивању нанотехнолошких уређаја који нису или не могу бити направљени.

У примењеној математици и теоријској физици, МД је постала део истраживања у домену динамичких система, ергодске теорије и статистичке механике. Принцип очувања енергије и ентропије потиче из термодинамике. Неке теорије за израчунавање конформационе ентропије, као што су основне компоненте анализе, односно техника која се користи за редукцију димензија скупа података ради њихове анализе, проистичу из информационе теорије. Математичке технике као што је трансформациони оператор, постају применљиве када МД посматрамо као Марковљев ланац. Постоји велики број математичара који користе симплектички интегратор, који представља нумеричку интеграциону шему, ради побољшања и ефикасности МД симулација. Симплектички интегратор спада у класу геометријских интегратора и предвиђен је за решавање специјалних диференцијалних једначина које су у вези са класичном механиком.

МД можемо посматрати као специјалан случај дискретне елементарне методе (ДЕМ), у коме честице имају сферни облик (нпр. величине Ван дер Валсовог радијуса).

Дизајн МД симулација одређен је могућим капацитетом компјутера. Величине симулације (n=број честица), временски корак и укупно време трајања симулације треба одабрати тако да се прорачун заврши у неком разумном временском периоду. Ипак, симулација треба да буде довољно дуга да би верно представила природни процес који проучавамо. Већина научних радова на тему динамике протеина и ДНК користе симулацију у временском интервалу од наносекунде (10-9 s) до микросекунде (10-6 s). Да би се одржале симулације, потребно је од неколико дана до неколико година. Паралелно алгоритми дозвољавају подацима да се распореде дуж CPU-a, на пример просторна декомпозиција у LAMMPS.

Током класичне МД симулације, најтежи CPU задатак је процена потенцијала (поље сила) у функцији унутрашњих координата честица. У оквиру процене енергије, најважнији је невезивни или нековалентни део. У Великој О нотацији (асимптотској нотацији), уобичајена молекулска стимулациона скала је О (n2), ако су сви електростатични парови и ван дер Валсове интеракције узете експлицитно у рачун. Овај поступак можемо скратити ако употребимо електростатичке методе као што су Евалдова сфера (O(nlog(n))) или сферну граничну технику (O(n)).

Још један фактор који утиче на укупно CPU време симулације је величина интеграције временског корака. То је временски период између мерења потенцијала. Временски корак мора бити довољно мали да би се избегле дискретизационе грешке. Типичан временски корак у класичној МД је реда величине фемтосекунде (10-15 s). Ова вредност може бити повећана коришћењем алгоритама као што су SHAKE i LINCS, који постављају најбрже атоме на место. Развијене су и методе променљивог временског корака.^[5][6][7]

Морамо направити селекцију између експлицитног и имплицитног раствора. Честице експлицитног раствора (као што су TIP3P i SPC/E модели воде)^[8] морају бити израчунате преко поља сила, док за имплицитне растворе користимо приступ "усредњеног поља" (mean-field, енг.). Утицај експлицитних раствора на време потребно за компјутерску обраду података-CPU време(CPU, central processing unit) може бити десетострук или већи. Величина честица и вискозност експлицитних раствора је суштинска за одређивање особина молекула раствора.

У симулацији молекула експлицитног раствора, експлицитни простор мора бити довољно велики да би се избегли гранични услови. Гранични услови се најчешће уклањају бирањем фиксних вредности на крајевима или употребом периодичних граничних услова, када се једна страна симулације везује за супротну страну, опонашајући фазу великог капацитета.

У микроканонском ансамблу не долази до промене броја честица (N), запремине (V) ни енергије (E). Овоме одговара адијабатска промена, без размене топлоте. Динамичка путања молекула одговара промени потенцијалне и кинетичке енергије, док укупна енергија остаје константна.

За систем од N честица, у Декартовом координатном систему једначине кретања су дате Њутновом механиком, као

${\displaystyle F(X)=-\nabla U(X)=M{\dot {V}}(t)}$

${\displaystyle V(t)={\dot {X}}(t)}$

Потенцијална енергију U(X) је функција парцијалне координате X, идентична потенцијалу у физици и пољу сила у хемији. Прва једначина проистиче из Њутнових закона; сила F која делује на сваку честицу у систему може се израчунати као негативни градијент потенцијала U(X).

За сваки временски корак, сваки положај X, и сваку брзину V можемо интегралити помоћу симплектичке методе, као што је Верлет. Промена X и V са временом назива се трајекторија. Узимајући у обзир почетне услове и брзину, можемо да израчунамо све будуће (или прошле) положаје и брзине.

Један од збуњујућих фактора у МД је значење температуре. Најчешће имамо искуства са макроскопском температуром која укључује велики број честица. Али, температура је статистичка величина. Тако да, ако имамо велики број атома статистичку температуру можемо добити из тренутне температуре, коју добијамо када кинетичку енергију изједначимо са nkBT/2, где је n-број степени слободе.

Температура расте упркос малом броју атома који користимо у МД симулацији. На пример, замислите симулацију раста бакарног филма на супстрату који садржи 500 атома и енергију од 100 eV. У реалном случају, енергија од 100 eV која потиче од исталожених атома, распоредила би се и пренела на велики број атома (>од 1010), и не би дошло до мерљиве промене температуре. Када има само 500 атома, супстрат у тренутку испарава услед разлагања. Нешто слично се дешава у биофизичкој симулацији. Температура система у NVE расте када макромолекул као што је протеин, подлеже егзотермној конформационој промени и везивању.

У канонском ансамблу количина супстанце N, запремина V i температура T су константни. Канонски ансамбл се понекад назива молекулска динамика на константној температури (ЦТМД). У NVT енергија ендотермног и егзотермног процеса размењује се са термостатом.

Да би додавање и одузимање енергије са граница МД система изгледало реалније, тј. апроксимирајући канонски ансамбл, условљава потребу низа термостатичких метода.

Популарне технике за контролисање температуре укључују Ноз-Хуверов термостат и Ланжеову динамику.

У NPT ансамблу количина супстанце N, притисак P и температура T су константни. При употреби термостата потребно је користити и баростат. Овакво стање највише одговара лабораторијским условима, отвореног суда на собној температури и притиску.

У симулацији биолошких мембрана, контрола изотропског притиска није задовољавајућа. За двослојне липиде контрола притиска се врши на константној површини мембране (NPAT) или на константном површинском напону, (NPγT).

Копија метода измене представља генерализован ансамбл. Његова главна сврха је била бављење спором динамиком неуређених спинских система. Изменом температура неинтерагујућих копија система при раду на више различитих температура, копија измене МД (РЕМД)^[9] на тај начин покушава да превазиђе вишеструке проблеме.

МД симулације захтевају дефиницију потенцијалне функције, или опис термова између којих ће честица у симулацији интераговати. У хемији и биологији под овим подразумевамо поље сила. Потенцијал може бити дефинисан на више начина; они који се обично користе у хемији, засновани су на молекулској механици, обухватају класичну теорију интеракције између две честице, које може да прикаже структурну и конформациону промену, али не и хемијску реакцију. Када треба дати детаљнији опис, користи се потенцијал дефинисан квантном механиком; поједине технике покушавају да створе хибрид класично-квантни потенцијал, где се већи део посматра класично а мањи део квантно механички, који обично подлеже хемијској трансформацији.

Већина поља сила у хемији су емпиријска, која се састоје од суме ковалентних веза које су у релацији са хемијским везама, угловима између њих и угловима између двоструких веза раздвојене трећом везом, и нековалентне везе које су повезане са ван дер Валсовим силама и наелектрисањем. Емпиријски потенцијали представљају квантно-механички ефекат, ограничени само за овакве функцијске апроксимације. Ови потенцијали садрже слободне параметре као што су наелектрисање атома, ван дер Валсови параметри као што су радијус атома, дужина везе; они су добијени електронским калкулацијама (квантно-хемијским симулацијама), или преко физичких особина као што су еластична константа, параметри решетке и спектроскопска мерења. Емпиријски потенцијали обично могу да представе поједину потенцијалну површ, и тако могу да прикажу процес раскидања хемијских веза, као и саме хемијске везе. Ово је директно последица коришћења Борн-Опенхајмерове апроксимације, или адијабатског формализма. Ако потенцијал одговара основном стању, не може се променити.

Потенцијалне функције које приказују невезивне интеракције, обично су "пар потенцијал", у коме се укупна енергија система рачуна као сума свих парцијалних чланова енергије сваког атома. Ове интеракције, пошто су нелокализоване и укључују најслабије интеракције између сваког пара честица у систему, обично представљају сметњу у брзини МД симулације. Код електростатичких интеракција, решавање Поасонове једначине за цео систем је углавном споро и уместо тога користе се нумеричке апроксимације, као што су променљиви гранични радијуси, реакциона поља алгоритама, сумација честичне Евалдове сфере, или новија Честица-Честица Честица решетка, која представља специјалан хибридни компјутерски алгоритам, уз коришћење N тачака у симулацији (P3M). Пример израчунатог потенцијала за пар честица је Ленард-Џоунсов потенцијал (познат ко 6-12 потенцијал) који користимо за израчунавање ван дер Валсових сила.

${\displaystyle U(r)=4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]}$

Други пример је Борнов (јонски) модел јонске решетке. Први израз у следећој једначини је Кулонов закон за јонски пар, а други израз је објашњен Паулијевим принципом и последњи израз је дисперзиона интеракција. У симулацију се укључује обично само диполарни терм (понекад квадрополни).

${\displaystyle U_{ij}(r_{ij})=\sum {\frac {z_{i}z_{j}}{4\pi \epsilon _{0}}}{\frac {1}{r_{ij}}}+\sum A_{l}exp{\frac {-r_{ij}}{p_{l}}}+\sum C_{l}r_{ij}^{-n_{j}}+...}$

Емпиријски потенцијали могу да се поделе на пар потенцијала и на потенцијал више честица. Код потенцијал који потиче од вишечестичног система, потенцијална енергија се не може израчунати као сума парова атома. На пример, Терсоф потенцијал, који се користи за симулацију силицијума и германијума, укључује суму група од три атома са угловима између атома као важним фактором потенцијала. Други пример је "Уско-Везивање Другог Момента Апроксимације" (Tight-Binding Second Moment Approximation, TBSMA, eng.) где електронску густину атома рачунамо из суме доприноса свих околних атома, а потенцијална енергија је тада функција те суме.

Семи-емпиријски потенцијал је користан у матричној репрезентацијама квантне механике. Ипак, вредности матричних елемената се добијају помоћу емпиријских формула које рачунају степен прелаза у одређеним атомским орбиталама. Рачунамо дијагоналу матрице да утврдимо попуњеност различитих атомских орбитала, и поново користимо емпиријске формуле да утврдимо енергијски допринос орбиталама.

Постоји пуно семи-емпиријских потенцијала, познатих као уско-везани потенцијали који се мењају, судећи на основу моделирања атома.

Ab - initio методе (први принципи) користе у потпуности квантно-механичку формулу за израчунавање потенцијалне енергије атома или молекула у систему. У поређењу са класичном потенцијалном функцијом која представљамо емпиријским функцијама, особине система представљамо рачунајући таласне функције за електроне који се крећу око језгра атома. Овај рачун се обично примењује локално, на пример за нуклеус у близини реакционе координате. Иако можемо да употребимо различите апроксимације, ове су базиране на теоријским претпоставкама, а не на емпиријским законима. Ab - initio метода обезбеђује пуно информација које не можемо добити емпиријским путем. Значајна предност ових метода је могућност проучавања реакција које укључују раскидање или стварање ковалентних веза, што одговара вишеструким електронским стањима. Класична молекуларна динамика не може да представи раскидање и формирање ковалентних веза. Ипак, последњих година технике као што су термодинамичка интеграција и дух честице, полако савладавају ове препреке. Успех је ипак ограничен.

Популаран пакет Ab - initio МД је Кар-Паринелов пакет (CPMD), заснован на теорији функционалне густине.

QM (квантно-механичке) методе су врло моћне, мада су врло скупе, док су ММ (класичне или молекулско механичке) методе брзе, али имају пар ограничења (захтевају екстензивну параметризацију, процена енергије није потпуно тачна, не могу се користити за симулацију реакција код којих се раскидају (или формирају) ковалентне везе, и не могу да обезбеде тачне детаље о хемијској околини). Нова класа метода обједињује предности QM (тачност) и ММ (брзину) метода. Ове методе су познате као хибрид QM и MM. Методологију за ове технике је представио А. Варшел са сарадницима. Последњих година појавило се неколико група које се баве овом методом, као што су: А. Варшал (Универзитет у Јужној Калифорнији), Кенет Мерз (Универзитет на Флориди), и други.

Најважнија предност хибрида QM/MM је брзина, док рад са класичном молекулском динамиком која непосредно користи скалу O (n²), где је N број атома у систему, представља главни проблем у изразу за електростатичку интеракцију (свака честица интерагује са осталим честицама у систему). Мада, коришћење граничног радијуса, повремене измене и најскорије промене у методи Евалдове честичне сфере, смањиле су скалу од O (n) до O (n²). Другим речима, ако симулирамо систем са дупло више атома, биће потребан од два до четири пута већи капацитет компјутера. Са друге стране, најједноставнији аб-инитио прорачуни су на скали O (n³), или као што је Хартри-Фок прорачун био на скали ~O (nn^2.7). Да би се ово превазишло, мали део система се третира квантно механички, а преостали део система класично.

У софистициранијој примени, QM/MM методе постоје да би објасниле сусцептибилност лаког језгра квантним ефектом (као што је водоник) и електронских стања, што дозвољава извођење водоникове таласне функције (слично електронској таласној функцији). Ова методологија је била успешна у објашњавању феномена, као што је тунел ефекат код водоника. Један пример где су ове методе дале нова сазнања је код прорачуна трансфера хидрида у ензиму јетре, алкохол-дехидрогеназе. У овом случају, тунел ефекат је важан за водоник да би се утврдила брзина реакције.^[10]

На другом крају скале су крупне честице и модел решетке. Уместо да експлицитно представљамо сваки атом система, користићемо псеудо атоме за представљање групе атома. МД симулације за велике системе захтевају велике компјутерске могућности тако да се не могу лако представити традиционалним методама у којима су приказани сви атоми. Слично, симулације процеса дужих од 1 µс су веома скупе, јер захтевају много међукорака. Овај проблем можемо решити користећи редуковане репрезентације, које се још називају и модели крупне-честице.

Примери метода крупне-честице (ЦГ) су променљива молекулска динамика (ЦГ-ДМД)^[11] и Го-модели.^[12] Понекад ЦГ користи пуно псеудо атома, које је као апроксимацију користила МД у симулацији биолошке мембране. Алифатични крај (реп) липида представљамо са неколико псеудоатома, тако што се 2-4 метиленске групе групишу у један псеудоатом.

Параметризација ЦГ модела мора да се уради емпиријски, усклађујући понашање модела са одговарајућим експерименталним подацима или симулацијом у којој учествују сви атоми. У идеалном случају, ови параметри и за енталпију и за ентропијски имплицитни допринос слободној енергији. Када се ЦГ изврши за више нивое, тачност динамичког описа може бити мање поуздан. Али многи ЦГ модели су успешно коришћени при испитивању широког опсега питања у структуралној биологији.

Примери примене ЦГ у биофизици:

• проучавања увијања протеина су често извођена коришћењем једног или неколико псеудоатома по једној амино киселини
• спирално увијање ДНА ланца је истражено коришћењем 1-3 псеудоатома по основном пару, или чак и на нижој резолуцији
• паковање хеликса ДНА у бактериофагу је истраживано моделима где 1 псеудоатом представља један окрет (око 10 основних парова) хеликса
• РНА структура рибозома и других великих система је представљена једним псеудоатомом по нуклеотиду

Најједноставнији облик ЦГ је „уједињен атом“ (понекад називан „продужен атом“) и примењиван у најранијим МД симулацијама протеина, липида и нуклеинских киселина. На пример, уместо да посматрамо сва четири атома метил групе CH3(или сва три атома метиленске групе CH2), један псеудоатом представља целу групу. Он мора бити добро окарактерисан тако да ван дер Валсове интеракције са другим групама буду зависне на одговарајућој удаљености. Слично важи за везе, углове и торзионе углове, где примењујемо псеудоатоме. Овој врсти репрезентације, где један карактеристичан елиминише све водоникове атоме, осим оних који су способни да граде водоничне везе. Пример овога је поље сила Чарм 19. Поларни водоникови атоми су обично задржани у моделу, јер добра представа водоничних веза захтева тачан опис усмерености и електростатичких интеракција између донора и акцептора група. Хидроксилна група, на пример може бити и донор и акцептор водоничне везе, и било би немогуће да се ово представи само једним OH псеудоатомом. Примећујемо да скоро половина атома у протеину или нуклеинској киселини представљају неполарни водоници, па употреба уједињених атома обезбеђује уштеду времена.

Молекулска динамика се користи у многим областима науке.

• Прва објављена симулација макромолекула (1977, величина:500 атома, време симулације: 9.2 ps, програм: Чарм процесор) протеин: говеђи панкреас трипсин инхибитор. Ово је једна од најбољих студија протеина по питању представе о увијању и кинетике протеина. Ова симулација објављена у магазину „Природа“(Nature, енг), прокрчила је пут ка разумевању кретања протеина, што представља суштину протеина а не само опис.^[13]

Следећа два биофизичка примера нису пресечне МД симулације. Они представљају смели подухват да се направи симулација веома великог система (комплетан вирус) који траје јако дуго(500μs):

• МД симулација СТМВ (2006, величина: 1.000.000 атома, време симулације: 50ns, програм: NAMD) овај вирус је мали биљни вирус који погоршава симптоме инфекције ТМВ вируса. МД симулације су коришћене за истраживање механизама вирусних група. Цео СТМВ се састоји од 60 идентичних копија једног протеина који прави спољашњи омотач вируса, и 1063 нуклеотида једног РНА генома. Један начин посматрања је да је спољашњи омотач нестабилан када унутар њега не постоји ланац РНА. За симулацију би било потребно 2006 компјутера и око 35 година да се заврши. То је урађено са много процесора који су упоредо радили, и били у међусобној вези.^[14]
• Симулација увијања домена протеина Вилина (Vilin headpiece) у програм: Folding@home 2006 (величина: 20.000 атома, време симулације: 500 μs). Ова симулација је рађена помоћу 200.000 личних рачунара широм света. Ови компјутери су имали инсталирани Folding@home програм, и било је координирано са Станфордског Универзитета. Кинетички параметри тог протеина били су истражени коришћењем независних, кратких трајекторија рађених помоћу CPU-a без континуалне временске комуникације. Једна коришћена техника је била Pfold анализа (Pfold је сервер за одређивање секвенце РНА)^[15] која даје информације о прелазним стањима структура и распореду конформација у току увијања. Свака трајекторија у Pfold прорачуну може бити релативно краћа, али потребне су многе независне трајекторије.^[16]

• ABINIT (ДФТ)
• AMBER (класична)
• Ascalaph (класична)
• CASTEP (ab-initio)
• CPMD (ДФТ)
• CHARMM (класична)
• CP2K (ДФТ)
• ESPResSo(класична, паралелна, продужена)
• Fireball (ab-initio)
• DL_POLY (класична)
• GROMOS (класична)
• GROMACS (класична)
• LAMMPS (класична, симулација паралелног посматрања са високом скалом просторне декомпозиције)
• NAMD (класична)
• ProtoMol (класична)
• PWscf
• SIESTA(ДФТ)
• VASP (ab-initio)
• TINKER (класична)