Молекулярная динамика ( МД ) – это метод компьютерного моделирования для анализа физических движений атомов и молекул . Атомам и молекулам разрешено взаимодействовать в течение фиксированного периода времени, что дает представление о динамической «эволюции» системы. В наиболее распространенной версии траектории атомов и молекул определяются путем численного решения уравнений движения Ньютона для системы взаимодействующих частиц, где силы между частицами и их потенциальные энергии часто рассчитываются с использованием межатомных потенциалов или молекулярных механических силовых полей . Метод применяется главным образом в химической физике , материаловедении и биофизике .
Поскольку молекулярные системы обычно состоят из огромного числа частиц, аналитически определить свойства таких сложных систем невозможно ; МД-моделирование позволяет обойти эту проблему, используя численные методы. Однако длительное МД-моделирование математически плохо обусловлено , генерируя совокупные ошибки численного интегрирования, которые можно минимизировать при правильном выборе алгоритмов и параметров, но не устранить.
Для систем, которые подчиняются эргодической гипотезе , эволюция одного моделирования молекулярной динамики может быть использована для определения макроскопических термодинамических свойств системы: средние значения по времени эргодической системы соответствуют средним значениям микроканонического ансамбля . MD также называют « статистической механикой в числах» и « видением Лапласа ньютоновской механики », предсказывающей будущее путем оживления сил природы [1] и позволяющей понять молекулярное движение в атомном масштабе.
Первоначально МД была разработана в начале 1950-х годов, после более ранних успехов моделирования Монте-Карло , которые сами по себе относятся к восемнадцатому веку, например, в задаче об игле Бюффона , но была популяризирована для статистической механики в Национальной лаборатории Лос-Аламоса Розенблутом и Метрополисом. в том, что сегодня известно как алгоритм Метрополиса – Гастингса . Интерес к временной эволюции систем N тел возник гораздо раньше, в семнадцатом веке, начиная с Ньютона, и продолжался в следующем столетии, в основном с акцентом на небесную механику и такие вопросы, как стабильность Солнечной системы. Многие из численных методов, используемых сегодня, были разработаны в этот период времени, до использования компьютеров; например, наиболее распространенный алгоритм интегрирования, используемый сегодня, алгоритм интегрирования Верле , был использован еще в 1791 году Жаном Батистом Жозефом Деламбром . Численные расчеты с использованием этих алгоритмов можно считать МД «вручную».
Еще в 1941 году интегрирование уравнений движения многих тел проводилось с помощью аналоговых компьютеров. Некоторые взялись за трудоемкую работу по моделированию движения атомов путем построения физических моделей, например, с использованием макроскопических сфер. Целью было расположить их таким образом, чтобы воспроизвести структуру жидкости и использовать ее для изучения ее поведения. Дж. Д. Бернал сказал в 1962 году: « ... я взял несколько резиновых шариков и склеил их стержнями разной длины от 2,75 до 4 дюймов. Прежде всего я старался сделать это как можно более непринужденно. , работая в своем собственном кабинете, меня прерывают примерно каждые пять минут и я не помню , что я делал до перерыва » .
После открытия микроскопических частиц и развития компьютеров интерес расширился за пределы испытательного полигона гравитационных систем к статистическим свойствам материи. В попытке понять происхождение необратимости Ферми предложил в 1953 году и опубликовал в 1955 году [3] использование MANIAC I , также находящегося в Национальной лаборатории Лос-Аламоса , для решения временной эволюции уравнений движения для многих система тела подчиняется нескольким вариантам силовых законов; сегодня эта плодотворная работа известна как проблема Ферми-Пасты-Улама-Цингу . Эволюция энергии исходной работы во времени показана на рисунке справа.
В 1957 году Алдер и Уэйнрайт [4] использовали компьютер IBM 704 для моделирования идеально упругих столкновений между твердыми сферами . [4] В 1960 году, возможно, в первой реалистичной модели материи, Гибсон и др. смоделировал радиационное повреждение твердой меди , используя отталкивающее взаимодействие типа Борна-Майера вместе с поверхностной силой сцепления. [5] В 1964 году Рахман [6] опубликовал моделирование жидкого аргона , в котором использовался потенциал Леннарда-Джонса ; расчеты свойств системы, таких как коэффициент самодиффузии , хорошо сравниваются с экспериментальными данными. [6] Сегодня потенциал Леннарда-Джонса по-прежнему является одним из наиболее часто используемых межмолекулярных потенциалов : [7] [8] Он используется для описания простых веществ (также известных как Леннард-Джонсий [9] [10] [11] ) для концептуальные и модельные исследования, а также в качестве строительного блока во многих силовых полях реальных веществ. [12] [13]
Впервые использованный в теоретической физике , метод МД вскоре приобрел популярность в материаловедении , а с 1970-х годов он также стал широко применяться в биохимии и биофизике . MD часто используется для уточнения трехмерных структур белков и других макромолекул на основе экспериментальных ограничений, полученных из рентгеновской кристаллографии или ЯМР-спектроскопии . В физике МД используется для изучения динамики явлений атомного уровня, которые невозможно наблюдать напрямую, таких как рост тонких пленок и ионная субплантация, а также для изучения физических свойств нанотехнологических устройств, которые еще не созданы или еще не могут быть созданы. . В биофизике и структурной биологии этот метод часто применяется для изучения движений макромолекул, таких как белки и нуклеиновые кислоты , что может быть полезно для интерпретации результатов определенных биофизических экспериментов и для моделирования взаимодействий с другими молекулами, например, при стыковке лигандов . В принципе, МД можно использовать для предсказания структуры белка ab initio путем моделирования сворачивания полипептидной цепи из случайного клубка .
Результаты МД-моделирования можно проверить путем сравнения с экспериментами по измерению молекулярной динамики, популярным методом которых является ЯМР-спектроскопия. Прогнозы структуры, полученные с помощью MD, можно проверить с помощью экспериментов в масштабах всего сообщества по критической оценке прогнозирования структуры белка ( CASP ), хотя исторически этот метод имел ограниченный успех в этой области. Майкл Левитт , получивший Нобелевскую премию частично за применение МД к белкам, писал в 1999 году, что участники CASP обычно не использовали этот метод из-за « ... центрального затруднения молекулярной механики, а именно того, что минимизация энергии или молекулярная динамика в целом приводит к модели , которая меньше похожа на экспериментальную структуру » . доработка, не доходя до практической полезности в этих областях; многие считают параметры силового поля ключевой областью для дальнейшего развития. [15] [16] [17]
Сообщалось о МД-моделировании для разработки фармакофоров и дизайна лекарств. [18] Например, Пинто и др. реализовали МД-моделирование комплексов Bcl-Xl для расчета средних положений критических аминокислот, участвующих в связывании лигандов. [19] С другой стороны, Карлсон и др. реализовали моделирование молекулярной динамики для идентификации соединений, которые дополняют рецептор, вызывая при этом минимальное нарушение конформации и гибкости активного сайта. Снимки белка через постоянные интервалы времени во время моделирования накладывались для идентификации консервативных областей связывания (консервативных как минимум в трех из одиннадцати кадров) для разработки фармакофора. Спиракис и др. полагались на рабочий процесс MD-моделирования, отпечатков пальцев для лигандов и белков (FLAP) и линейного дискриминантного анализа для определения лучших конформаций лиганд-белок, которые будут действовать в качестве шаблонов фармакофоров, на основе ретроспективного ROC-анализа полученных фармакофоров. Пытаясь улучшить моделирование открытия лекарств на основе структуры, учитывая необходимость во многих смоделированных соединениях, Hatmal et al. предложили сочетание МД-моделирования и анализа межмолекулярных контактов лиганд-рецептор, чтобы отличить критические межмолекулярные контакты (связывающие взаимодействия) от избыточных в одном лиганд-белковом комплексе. Критические контакты затем могут быть преобразованы в модели фармакофора, которые можно использовать для виртуального скрининга. [20]
Важным фактором являются внутримолекулярные водородные связи , [21] которые не включены явно в современные силовые поля, а описываются как кулоновские взаимодействия атомных точечных зарядов. [ нужна цитата ] Это грубое приближение, поскольку водородные связи имеют частично квантово-механическую и химическую природу. Кроме того, электростатические взаимодействия обычно рассчитываются с использованием диэлектрической проницаемости вакуума , хотя окружающий водный раствор имеет гораздо более высокую диэлектрическую проницаемость . Использование макроскопической диэлектрической проницаемости на малых межатомных расстояниях сомнительно. Наконец, взаимодействия Ван-дер-Ваальса в МД обычно описываются потенциалами Леннарда-Джонса, основанными на теории Фрица Лондона , применимой только в вакууме. [ нужна цитация ] Однако все виды сил Ван-дер-Ваальса в конечном итоге имеют электростатическое происхождение и поэтому зависят от диэлектрических свойств окружающей среды . [22] Прямое измерение сил притяжения между различными материалами (как константа Гамакера ) показывает, что «взаимодействие между углеводородами в воде составляет около 10% от взаимодействия в вакууме». [22] Зависимость сил Ван-дер-Ваальса от окружающей среды не учитывается в стандартном моделировании, но может быть учтена путем разработки поляризуемых силовых полей.
При разработке моделирования молекулярной динамики следует учитывать доступную вычислительную мощность. Размер моделирования ( n = количество частиц), временной шаг и общая продолжительность времени должны быть выбраны так, чтобы расчет мог завершиться в течение разумного периода времени. Однако моделирование должно быть достаточно продолжительным, чтобы соответствовать временным масштабам изучаемых природных процессов. Чтобы сделать статистически обоснованные выводы из моделирования, моделируемый промежуток времени должен соответствовать кинетике естественного процесса. В противном случае это аналогично выводу о том, как человек ходит, глядя только на один шаг. В большинстве научных публикаций о динамике белков и ДНК [23] [24] используются данные моделирования, охватывающие периоды от наносекунд (10–9 с ) до микросекунд (10–6 с ). Для получения таких симуляций требуется от нескольких ЦП-дней до ЦП-лет. Параллельные алгоритмы позволяют распределять нагрузку между процессорами; примером является алгоритм пространственной или силовой декомпозиции. [25]
Во время классического МД-моделирования наиболее ресурсоемкой задачей ЦП является оценка потенциала как функции внутренних координат частиц. В рамках этой энергетической оценки самой дорогой является несвязанная или нековалентная часть. В обозначении Big O общие симуляции молекулярной динамики масштабируются в зависимости от того, необходимо ли явно учитывать все парные электростатические и ван-дер-ваальсовые взаимодействия . Эти вычислительные затраты можно уменьшить, используя методы электростатики, такие как суммирование сетки частиц Эвальда ( ), частица-частица-частица-сетка ( P3M ) или хорошие методы сферического отсечения ( ). [ нужна цитата ]
Еще одним фактором, влияющим на общее время процессора, необходимое для моделирования, является размер временного шага интегрирования. Это промежуток времени между оценками потенциала. Временной шаг должен быть выбран достаточно малым, чтобы избежать ошибок дискретизации (т. е. меньшим, чем период, относящийся к самой высокой частоте колебаний в системе). Типичные временные шаги для классической МД составляют порядка 1 фемтосекунды (10–15 с ). Это значение можно увеличить с помощью таких алгоритмов, как алгоритм ограничения SHAKE , который фиксирует колебания самых быстрых атомов (например, водородов) на месте. Также были разработаны методы с несколькими временными шкалами, которые позволяют увеличить время между обновлениями более медленных сил дальнего действия. [26] [27] [28]
Для моделирования молекул в растворителе необходимо сделать выбор между явным и неявным растворителем . Частицы явного растворителя (такие как модели воды TIP3P , SPC/E и SPC-f ) должны рассчитываться с помощью силового поля, в то время как неявные растворители используют подход среднего поля. Использование явного растворителя требует больших вычислительных затрат и требует включения в моделирование примерно в десять раз больше частиц. Но зернистость и вязкость явного растворителя необходимы для воспроизведения определенных свойств молекул растворенного вещества. Это особенно важно для воспроизведения химической кинетики .
Во всех видах моделирования молекулярной динамики размер окна моделирования должен быть достаточно большим, чтобы избежать артефактов граничных условий . Граничные условия часто обрабатываются путем выбора фиксированных значений на краях (что может вызвать артефакты) или использования периодических граничных условий , при которых одна сторона моделирования возвращается к противоположной стороне, имитируя объемную фазу (что также может вызывать артефакты). .
В микроканоническом ансамбле система изолирована от изменений молей (N), объема (V) и энергии (E). Это соответствует адиабатическому процессу без теплообмена. Траекторию микроканонической молекулярной динамики можно рассматривать как обмен потенциальной и кинетической энергией с сохранением полной энергии. Для системы N частиц с координатами и скоростями следующую пару дифференциальных уравнений первого порядка можно записать в обозначениях Ньютона как
Потенциальная энергия системы является функцией координат частицы . Его называют просто потенциалом в физике или силовым полем в химии. Первое уравнение происходит из законов движения Ньютона ; силу , действующую на каждую частицу в системе, можно рассчитать как отрицательный градиент .
Для каждого временного шага положение и скорость каждой частицы можно проинтегрировать с помощью метода симплектического интегратора , такого как интегрирование Верле . Эволюция во времени и называется траекторией. Учитывая начальные положения (например, из теоретических знаний) и скорости (например, рандомизированные гауссовы), мы можем вычислить все будущие (или прошлые) положения и скорости.
Одним из частых источников путаницы является значение температуры в MD. Обычно у нас есть опыт работы с макроскопическими температурами, в которых участвует огромное количество частиц, но температура — это статистическая величина. Если число атомов достаточно большое, статистическую температуру можно оценить по мгновенной температуре , которая находится путем приравнивания кинетической энергии системы к nk BT /2, где n — число степеней свободы системы.
Явление, связанное с температурой, возникает из-за небольшого количества атомов, которые используются в МД-моделировании. Например, рассмотрим моделирование роста медной пленки, начиная с подложки, содержащей 500 атомов и энергии осаждения 100 эВ. В реальном мире энергия 100 эВ от осажденного атома быстро переносится и распределяется между большим количеством атомов ( или более) без большого изменения температуры. Однако когда атомов всего 500, подложка практически сразу испаряется при осаждении. Нечто подобное происходит и в биофизическом моделировании. Температура системы при NVE естественным образом повышается, когда макромолекулы, такие как белки, подвергаются экзотермическим конформационным изменениям и связыванию.
В каноническом ансамбле сохраняются количество вещества (N), объём (V) и температура (Т). Его также иногда называют молекулярной динамикой при постоянной температуре (CTMD). В НВТ происходит обмен энергией эндо- и экзотермических процессов с термостатом.
Доступны различные алгоритмы термостата для более или менее реалистичного добавления и удаления энергии за границы МД-моделирования, аппроксимируя канонический ансамбль . Популярные методы контроля температуры включают изменение масштаба скорости, термостат Нозе-Гувера , цепи Нозе-Гувера, термостат Берендсена , термостат Андерсена и динамику Ланжевена . Термостат Берендсена может вызвать эффект летающего кубика льда , который приводит к нефизическим перемещениям и вращениям моделируемой системы.
С помощью этих алгоритмов получить каноническое ансамблевое распределение конформаций и скоростей непросто . Как это зависит от размера системы, выбора термостата, параметров термостата, временного шага и интегратора, является предметом многих статей в этой области.
В изотермо-изобарном ансамбле сохраняются количество вещества (N), давление (Р) и температура (Т). Помимо термостата необходим баростат. Он наиболее точно соответствует лабораторным условиям с колбой, открытой при температуре и давлении окружающей среды.
При моделировании биологических мембран изотропное управление давлением нецелесообразно. Для липидных бислоев контроль давления происходит при постоянной площади мембраны (NPAT) или постоянном поверхностном натяжении «гамма» (NPγT).
Метод обмена репликами представляет собой обобщенный ансамбль. Первоначально он был создан для изучения медленной динамики неупорядоченных спиновых систем. Его еще называют параллельным отпуском. Формулировка MD с обменом репликами (REMD) [29] пытается преодолеть проблему нескольких минимумов путем обмена температур невзаимодействующих реплик системы, работающей при нескольких температурах.
Моделирование молекулярной динамики требует определения потенциальной функции или описания условий, с помощью которых частицы в моделировании будут взаимодействовать. В химии и биологии это обычно называют силовым полем , а в физике материалов — межатомным потенциалом . Потенциалы могут быть определены на многих уровнях физической точности; те, которые наиболее часто используются в химии, основаны на молекулярной механике и воплощают классическую механическую трактовку взаимодействий частица-частица, которая может воспроизводить структурные и конформационные изменения , но обычно не может воспроизводить химические реакции .
Сведение от полностью квантового описания к классическому потенциалу предполагает два основных приближения. Первое — это приближение Борна-Оппенгеймера , которое утверждает, что динамика электронов настолько быстра, что можно считать, что они мгновенно реагируют на движение своих ядер. В результате их можно рассматривать отдельно. Второй рассматривает ядра, которые намного тяжелее электронов, как точечные частицы, подчиняющиеся классической ньютоновской динамике. В классической молекулярной динамике влияние электронов аппроксимируется одной поверхностью потенциальной энергии, обычно представляющей основное состояние.
Когда необходимы более тонкие уровни детализации, используются потенциалы, основанные на квантовой механике ; некоторые методы пытаются создать гибридные классические/квантовые потенциалы, где основная часть системы рассматривается классически, но небольшая область рассматривается как квантовая система, обычно подвергающаяся химическому превращению.
Эмпирические потенциалы, используемые в химии, часто называют силовыми полями , а в физике материалов — межатомными потенциалами .
Большинство силовых полей в химии являются эмпирическими и состоят из суммы связанных сил, связанных с химическими связями , валентными углами и двугранниками связей , а также несвязанных сил, связанных с силами Ван-дер-Ваальса и электростатическим зарядом . [30] Эмпирические потенциалы представляют квантово-механические эффекты ограниченным образом посредством специальных функциональных приближений. Эти потенциалы содержат свободные параметры, такие как атомный заряд , параметры Ван-дер-Ваальса, отражающие оценки атомного радиуса , а также равновесную длину связи , угол и двугранник; они получены путем сопоставления подробных электронных расчетов (квантово-химическое моделирование) или экспериментальных физических свойств, таких как упругие константы , параметры решетки и спектроскопические измерения.
Из-за нелокальной природы несвязанных взаимодействий они включают, по крайней мере, слабые взаимодействия между всеми частицами в системе. Его расчет обычно является узким местом в скорости моделирования МД. Чтобы снизить вычислительные затраты, силовые поля используют численные аппроксимации, такие как смещенные радиусы отсечки, алгоритмы поля реакции , суммирование сетки частиц Эвальда или более новая сетка частица-частица-частица-сетка ( P3M ).
В химических силовых полях обычно используются заранее заданные схемы связей (исключением является динамика ab initio ), и поэтому они не могут явно моделировать процесс разрыва химических связей и реакции. С другой стороны, многие из потенциалов, используемых в физике, например, основанные на формализме порядка связей, могут описывать несколько различных координаций системы и разрыва связей. [31] [32] Примеры таких потенциалов включают потенциал Бреннера [33] для углеводородов и его дальнейшее развитие для систем C-Si-H [34] и COH [35] . Потенциал ReaxFF [36] можно рассматривать как полностью реактивный гибрид между потенциалами порядка связей и химическими силовыми полями .
Потенциальные функции, представляющие несвязанную энергию, формулируются как сумма по взаимодействиям между частицами системы. Самый простой выбор, используемый во многих популярных силовых полях , - это «парный потенциал», в котором полная потенциальная энергия может быть рассчитана из суммы энергетических вкладов между парами атомов. Поэтому эти силовые поля также называют «аддитивными силовыми полями». Примером такого парного потенциала является несвязанный потенциал Леннарда-Джонса (также называемый потенциалом 6–12), используемый для расчета сил Ван-дер-Ваальса.
Другой пример — борновская (ионная) модель ионной решетки. Первый член в следующем уравнении представляет собой закон Кулона для пары ионов, второй член — это короткодействующее отталкивание, объясняемое принципом Паули, а последний член — это член дисперсионного взаимодействия. Обычно моделирование включает только диполярный член, хотя иногда включается и квадрупольный член. [37] [38] Когда n l = 6, этот потенциал также называют потенциалом Кулона – Букингема .
В потенциалах многих тел потенциальная энергия включает в себя эффекты трех или более частиц, взаимодействующих друг с другом. [39] В симуляциях с парными потенциалами глобальные взаимодействия в системе также существуют, но они происходят только через парные члены. В потенциалах многих тел потенциальная энергия не может быть найдена путем суммы по парам атомов, поскольку эти взаимодействия рассчитываются явно как комбинация членов более высокого порядка. С статистической точки зрения зависимость между переменными вообще не может быть выражена с использованием только попарных произведений степеней свободы. Например, потенциал Терсоффа [40] , который первоначально использовался для моделирования углерода , кремния и германия и с тех пор использовался для широкого круга других материалов, включает в себя сумму по группам из трех атомов с углами между Атомы являются важным фактором потенциала. Другими примерами являются метод внедренного атома (EAM), [41] EDIP, [39] и потенциалы приближения второго момента сильной связи (TBSMA), [42] где плотность электронных состояний в области атома равна рассчитывается из суммы вкладов окружающих атомов, и тогда вклад потенциальной энергии является функцией этой суммы.
Полуэмпирические потенциалы используют матричное представление квантовой механики. Однако значения матричных элементов находятся с помощью эмпирических формул, оценивающих степень перекрытия конкретных атомных орбиталей. Затем матрица диагонализируется для определения занятости различных атомных орбиталей, а эмпирические формулы снова используются для определения энергетических вкладов орбиталей.
Существует большое разнообразие полуэмпирических потенциалов, называемых потенциалами сильной связи , которые варьируются в зависимости от моделируемых атомов.
Большинство классических силовых полей неявно включают эффект поляризуемости , например, путем увеличения парциальных зарядов, полученных в результате квантово-химических расчетов. Эти частичные заряды стационарны по отношению к массе атома. Но моделирование молекулярной динамики может явно моделировать поляризуемость с помощью индуцированных диполей с помощью различных методов, таких как частицы Друде или флуктуирующие заряды. Это позволяет осуществлять динамическое перераспределение заряда между атомами, которое реагирует на локальную химическую среду.
В течение многих лет поляризуемые модели МД рекламировались как следующее поколение. Для однородных жидкостей, таких как вода, повышенная точность была достигнута за счет учета поляризуемости. [43] [44] [45] Некоторые многообещающие результаты были также достигнуты для белков. [46] [47] Однако до сих пор неясно, как лучше всего аппроксимировать поляризуемость в моделировании. [ нужна цитата ] Этот вопрос становится более важным, когда частица испытывает различные условия во время своей траектории моделирования, например, транслокация лекарства через клеточную мембрану. [48]
В классической молекулярной динамике в силовом поле представлена одна поверхность потенциальной энергии (обычно основное состояние). Это следствие приближения Борна-Оппенгеймера . В возбужденных состояниях, химических реакциях или когда необходимо более точное представление, поведение электрона может быть получено из первых принципов с использованием квантово-механического метода, такого как теория функционала плотности . Это называется Ab Initio Molecular Dynamics (AIMD). Из-за стоимости обработки электронных степеней свободы вычислительная нагрузка этого моделирования намного выше, чем в классической молекулярной динамике. По этой причине AIMD обычно ограничивается меньшими системами и более коротким временем.
Квантово-механические и химические методы ab initio могут использоваться для расчета потенциальной энергии системы на лету, что необходимо для конформаций на траектории. Этот расчет обычно производится в непосредственной близости от координаты реакции . Хотя могут использоваться различные приближения, они основаны на теоретических соображениях, а не на эмпирических данных. Расчеты ab initio дают огромное количество информации, недоступной с помощью эмпирических методов, таких как плотность электронных состояний или другие электронные свойства. Существенным преимуществом использования методов ab initio является возможность изучения реакций, связанных с разрывом или образованием ковалентных связей, соответствующих множественным электронным состояниям. Более того, методы ab initio также позволяют восстанавливать эффекты, выходящие за рамки приближения Борна-Оппенгеймера, используя такие подходы, как смешанная квантово-классическая динамика .
КМ (квантово-механические) методы очень эффективны. Однако они являются дорогостоящими в вычислительном отношении, в то время как методы ММ (классической или молекулярной механики) быстры, но имеют несколько ограничений (требуют обширной параметризации; полученные оценки энергии не очень точны; не могут использоваться для моделирования реакций, в которых ковалентные связи разрываются/образуются). и ограничены в своих возможностях предоставлять точную информацию о химической среде). Появился новый класс методов, сочетающий в себе преимущества расчетов QM (точность) и MM (скорость). Эти методы называются смешанными или гибридными методами квантовой механики и молекулярной механики (гибридными КМ/ММ). [49]
Самым важным преимуществом гибридного метода QM/MM является скорость. Стоимость проведения классической молекулярной динамики (ММ) в наиболее простом случае масштабируется O(n 2 ), где n — количество атомов в системе. В основном это связано с термином электростатических взаимодействий (каждая частица взаимодействует с каждой другой частицей). Однако использование радиуса отсечки, периодические обновления списка пар и, в последнее время, вариации метода сетки частиц Эвальда (PME) снизили это значение до уровня от O(n) до O(n 2 ) . Другими словами, если смоделировать систему с вдвое большим количеством атомов, то потребуется в два-четыре раза больше вычислительной мощности. С другой стороны, простейшие расчеты ab initio обычно масштабируются O(n 3 ) или хуже ( было предложено ограничить масштаб вычислений Хартри-Фока для масштабирования ~ O (n 2,7 )). Чтобы преодолеть ограничение, небольшая часть системы обрабатывается квантово-механически (обычно активный центр фермента), а оставшаяся система обрабатывается классически.
В более сложных реализациях методы КМ/ММ существуют для обработки как легких ядер, чувствительных к квантовым эффектам (таких как водороды), так и электронных состояний. Это позволяет генерировать волновые функции водорода (аналогичные волновым функциям электронов). Эта методология оказалась полезной при исследовании таких явлений, как туннелирование водорода. Одним из примеров, когда методы КМ/ММ привели к новым открытиям, является расчет переноса гидрида в ферменте алкогольдегидрогеназы печени . В этом случае квантовое туннелирование важно для водорода, поскольку оно определяет скорость реакции. [50]
На другом конце шкалы детализации находятся крупнозернистые и решетчатые модели. Вместо явного представления каждого атома системы для обозначения групп атомов используются «псевдоатомы». МД-моделирование в очень больших системах может потребовать таких больших компьютерных ресурсов, что их невозможно легко изучить традиционными полноатомными методами. Точно так же моделирование процессов в длительных временных масштабах (более 1 микросекунды) является непомерно дорогим, поскольку требует очень большого количества временных шагов. В этих случаях иногда можно решить проблему, используя сокращенные представления, которые также называются крупнозернистыми моделями . [51]
Примерами методов грубого измельчения (CG) являются разрывная молекулярная динамика (CG-DMD) [52] [53] и Go-модели. [54] Иногда выполняется грубое измельчение с использованием более крупных псевдоатомов. Такие приближения объединенного атома использовались при МД-моделировании биологических мембран. Реализация такого подхода в системах, где электрические свойства представляют интерес, может быть сложной задачей из-за сложности использования правильного распределения заряда на псевдоатомах. [55] Алифатические хвосты липидов представлены несколькими псевдоатомами, собирая от 2 до 4 метиленовых групп в каждом псевдоатоме.
Параметризация этих очень грубых моделей должна осуществляться эмпирически, путем сопоставления поведения модели с соответствующими экспериментальными данными или полноатомным моделированием. В идеале эти параметры должны неявным образом учитывать как энтальпийный , так и энтропийный вклад в свободную энергию. [56] Когда грубая настройка выполняется на более высоких уровнях, точность динамического описания может быть менее надежной. Но очень грубые модели успешно использовались для изучения широкого круга вопросов структурной биологии, организации жидких кристаллов и полимерных стекол.
Примеры применения крупнозернистости:
Самая простая форма грубого анализа — это объединенный атом (иногда называемый расширенным атомом ), который использовался в большинстве ранних моделей МД белков, липидов и нуклеиновых кислот. Например, вместо того, чтобы явно рассматривать все четыре атома метильной группы CH 3 (или все три атома метиленовой группы CH 2 ), представляют всю группу с одним псевдоатомом. Разумеется, она должна быть правильно параметризована, чтобы ее взаимодействия Ван-дер-Ваальса с другими группами имели правильную зависимость от расстояния. Аналогичные соображения применимы и к связям, углам и кручениям, в которых участвует псевдоатом. В этом виде представления объединенного атома обычно исключаются все явные атомы водорода, за исключением тех, которые способны участвовать в водородных связях ( полярные водороды ). Примером этого является силовое поле CHARMM 19.
Полярные водороды обычно сохраняются в модели, поскольку правильное рассмотрение водородных связей требует достаточно точного описания направленности и электростатических взаимодействий между донорными и акцепторными группами. Гидроксильная группа, например, может быть как донором, так и акцептором водородной связи, и обработать ее одним псевдоатомом ОН было бы невозможно. Около половины атомов белка или нуклеиновой кислоты представляют собой неполярные атомы водорода, поэтому использование объединенных атомов может обеспечить существенную экономию компьютерного времени.
Во многих симуляциях системы растворенное вещество-растворитель основное внимание уделяется поведению растворенного вещества, при этом поведение растворителя мало интересует, особенно в тех молекулах растворителя, которые находятся в областях, далеких от молекулы растворенного вещества. [57] Растворители могут влиять на динамическое поведение растворенных веществ посредством случайных столкновений и создания сопротивления трения при движении растворенного вещества через растворитель. Использование непрямоугольных периодических граничных условий, стохастических границ и оболочек растворителя может помочь уменьшить количество требуемых молекул растворителя и позволить вместо этого тратить большую часть вычислительного времени на моделирование растворенного вещества. Также возможно включить эффекты растворителя без необходимости присутствия каких-либо явных молекул растворителя. Одним из примеров этого подхода является использование потенциальной средней силы (PMF), которая описывает, как изменяется свободная энергия при изменении конкретной координаты. Изменение свободной энергии, описываемое PMF, содержит усредненные эффекты растворителя.
Без учета эффектов растворителя моделирование макромолекул (таких как белки) может привести к нереалистичному поведению, и даже небольшие молекулы могут принять более компактные конформации из-за благоприятных сил Ван-дер-Ваальса и электростатических взаимодействий, которые будут ослабляться в присутствии растворителя. [58]
Дальнодействующее взаимодействие — это взаимодействие, при котором пространственное взаимодействие спадает не быстрее, чем размерность системы. Примеры включают взаимодействия заряд-заряд между ионами и диполь-дипольные взаимодействия между молекулами. Моделирование этих сил представляет собой довольно сложную задачу, поскольку они значительны на расстоянии, которое может превышать половину длины ящика при моделировании многих тысяч частиц. Хотя одним из решений могло бы стать значительное увеличение длины блока, этот подход грубой силы не идеален, поскольку моделирование стало бы очень дорогостоящим в вычислительном отношении. О сферическом усечении потенциала также не может быть и речи, поскольку может наблюдаться нереалистичное поведение, когда расстояние близко к предельному расстоянию. [59]
Моделирование управляемой молекулярной динамики (SMD), или моделирование силового зонда, применяет силы к белку, чтобы манипулировать его структурой, перемещая его по желаемым степеням свободы. Эти эксперименты можно использовать для выявления структурных изменений в белке на атомном уровне. SMD часто используется для моделирования таких событий, как механическое разворачивание или растяжение. [60]
Существует два типичных протокола SMD: один, в котором скорость вытягивания поддерживается постоянной, и другой, в котором приложенная сила постоянна. Обычно часть изучаемой системы (например, атом белка) удерживается гармоническим потенциалом. Затем к конкретным атомам прикладывают силы либо с постоянной скоростью, либо с постоянной силой. Зонтичная выборка используется для перемещения системы вдоль желаемой координаты реакции путем изменения, например, сил, расстояний и углов, которыми манипулируют при моделировании. Посредством зонтичного отбора проб все конфигурации системы – как высокоэнергетические, так и низкоэнергетические – подвергаются адекватной выборке. Затем изменение свободной энергии каждой конфигурации можно рассчитать как потенциал средней силы . [61] Популярным методом расчета PMF является метод анализа взвешенных гистограмм (WHAM), который анализирует серию моделей зонтичной выборки. [62] [63]
Многие важные применения SMD находятся в области открытия лекарств и биомолекулярных наук. Например, SMD использовался для исследования стабильности протофибрилл болезни Альцгеймера [64] для изучения взаимодействия белок-лиганд в циклин-зависимой киназе 5 [65] и даже для демонстрации влияния электрического поля на тромбин (белок) и аптамер (нуклеотид). комплекс [66] среди многих других интересных исследований.
Молекулярная динамика используется во многих областях науки.
Следующие биофизические примеры иллюстрируют заметные усилия по созданию моделирования систем очень большого размера (полный вирус) или очень длительного времени моделирования (до 1,112 миллисекунды):
Еще одно важное применение метода МД заключается в его способности трехмерной характеристики и анализа микроструктурной эволюции на атомном уровне.
Молекулярное моделирование на графическом процессоре — это метод использования графического процессора (GPU) для молекулярного моделирования. [79]
В 2007 году компания NVIDIA представила видеокарты, которые можно было использовать не только для отображения графики, но и для научных расчетов. Эти карты включают в себя множество арифметических единиц (по состоянию на 2016 год [update]— до 3584 в Tesla P100), работающих параллельно. Задолго до этого события вычислительные мощности видеокарт использовались исключительно для ускорения графических вычислений. Новым было то, что NVIDIA позволила разрабатывать параллельные программы в интерфейсе прикладного программирования (API) высокого уровня под названием CUDA . Эта технология существенно упростила программирование, позволив писать программы на C / C++ . Совсем недавно OpenCL позволил обеспечить кроссплатформенное ускорение графического процессора.