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