Молекулярная динамика: Компьютерное моделирование движения атомов и молекул.

Молекулярная динамика: компьютерное моделирование движения атомов и молекул для изучения химических свойств веществ. Применение в химии, физике и биологии.

25 мин чтения

100%

Источник: Molecular dynamics — Wikipedia

Сравнивайте с английским: нажмите на абзац — оригинал откроется в окне. Кнопка EN под абзацем показывает его прямо в тексте.

Содержание

Введение

Компьютерное моделирование для выявления и понимания химических свойств. Молекулярная динамика (МД) – это метод компьютерного моделирования, предназначенный для анализа физических перемещений атомов и молекул. Атомам и молекулам разрешается взаимодействовать в течение фиксированного периода времени, что позволяет наблюдать динамическую "эволюцию" системы. В наиболее распространенной реализации траектории атомов и молекул определяются численным решением уравнений движения Ньютона для системы взаимодействующих частиц, где силы между частицами и их потенциальные энергии часто рассчитываются с использованием межатомных потенциалов или молекулярно-механических силовых полей. Метод находит применение главным образом в химической физике, материаловедении и биофизике. Поскольку молекулярные системы обычно состоят из огромного числа частиц, аналитическое определение свойств таких сложных систем невозможно; МД-моделирование обходит эту проблему, используя численные методы. Однако длительные МД-моделирования математически неустойчивы, что приводит к накоплению ошибок при численном интегрировании, которые можно минимизировать, правильно выбирая алгоритмы и параметры, но не устранить полностью. Для систем, удовлетворяющих эргодической гипотезе, эволюция одного МД-моделирования может быть использована для определения макроскопических термодинамических свойств системы: временные средние эргодической системы соответствуют средним по микроканоническому ансамблю. МД также называют "статистической механикой, основанной на численном моделировании" и "ньютоновской механикой в представлении Лапласа", позволяющей предсказывать будущее, анимируя силы природы и давая представление о молекулярном движении на атомном уровне.

История

MD был первоначально разработан в начале 1950-х годов, после более ранних успехов с методами Монте-Карло, которые, в свою очередь, восходят к XVIII веку, например, к задаче Бюффона об игле, но получил широкое распространение в статистической механике в Лос-Аламосской национальной лаборатории благодаря работам Маршалла Розенблюта и Николаса Метрополиса, которые разработали алгоритм, известный сегодня как алгоритм Метрополиса — Хэстингса. Интерес к временной эволюции систем из N частиц возник гораздо раньше, в XVII веке, с работ Исаака Ньютона, и продолжился в следующем столетии, в основном с акцентом на небесной механике и вопросах, таких как стабильность Солнечной системы. Многие из численных методов, используемых сегодня, были разработаны в этот период, предшествующий появлению компьютеров; например, наиболее распространенный алгоритм интегрирования, используемый сегодня — алгоритм интегрирования Верле, был применен еще в 1791 году Жаном Батистом Жозефом Деламбром. Численные расчеты с использованием этих алгоритмов можно рассматривать как выполнение МД «вручную». Уже в 1941 году интегрирование уравнений движения для систем из многих частиц осуществлялось с помощью аналоговых компьютеров. Некоторые ученые взялись за трудоемкую задачу моделирования движения атомов, создавая физические модели, например, с использованием макроскопических сфер. Целью было расположить их таким образом, чтобы воспроизвести структуру жидкости и изучить ее поведение. Дж. Д. Берналь описал этот процесс в 1962 году, написав: «Я взял несколько резиновых шаров и склеил их стержнями различной длины, от 2,75 до 4 дюймов. Я старался делать это как можно более непринужденно, работая в своем кабинете, прерываясь примерно каждые пять минут и не помня, что я делал до перерыва». После открытия микроскопических частиц и развития компьютеров интерес расширился за пределы гравитационных систем и охватил статистические свойства материи. В попытке понять происхождение необратимости Энрико Ферми предложил в 1953 году и опубликовал в 1955 году использовать ранний компьютер MANIAC I, также в Лос-Аламосской национальной лаборатории, для решения задачи временной эволюции уравнений движения для системы из многих частиц, подверженной различным законам силы. Сегодня эта основополагающая работа известна как проблема Ферми — Паста — Улама — Цингоу. Временная эволюция энергии, полученная в оригинальной работе, показана на рисунке справа. В 1957 году Берни Алдер и Томас Уэйнрайт использовали компьютер IBM 704 для моделирования абсолютно упругих столкновений между твердыми сферами. В 1960 году, возможно, в первой реалистичной симуляции материи, Дж. Б. Гибсон и др. смоделировали радиационное повреждение твердой меди, используя отталкивающее взаимодействие типа Борна — Майера в сочетании с когезионной поверхностной силой. В 1964 году Анесур Рахман опубликовал результаты моделирования жидкого аргона, в котором использовался потенциал Леннарда — Джонса; расчеты свойств системы, такие как коэффициент самодиффузии, хорошо согласовались с экспериментальными данными. Сегодня потенциал Леннарда — Джонса по-прежнему является одним из наиболее часто используемых межмолекулярных потенциалов. Он используется для описания простых веществ (также известных как «Леннарда-Джонсиум») для концептуальных и модельных исследований, а также в качестве строительного блока во многих силовых полях для реальных веществ.

MD was originally developed in the early 1950s, following earlier successes with Monte Carlo simulationswhich themselves date back to the eighteenth century, in the Buffon's needle problem for examplebut was popularized for statistical mechanics at Los Alamos National Laboratory by Marshall Rosenbluth and Nicholas Metropolis in what is known today as the Metropolis–Hastings algorithm. Interest in the time evolution of N body systems dates much earlier to the seventeenth century, beginning with Isaac Newton, and continued into the following century largely with a focus on celestial mechanics and issues such as the stability of the solar system. Many of the numerical methods used today were developed during this time period, which predates the use of computers; for example, the most common integration algorithm used today, the Verlet integration algorithm, was used as early as 1791 by Jean Baptiste Joseph Delambre. Numerical calculations with these algorithms can be considered to be MD done "by hand". As early as 1941, integration of the many body equations of motion was carried out with analog computers. Some undertook the labor intensive work of modeling atomic motion by constructing physical models, e. g., using macroscopic spheres. The aim was to arrange them in such a way as to replicate the structure of a liquid and use this to examine its behavior. J. D. Bernal describes this process in 1962, writing: I took a number of rubber balls and stuck them together with rods of a selection of different lengths ranging from 2.75 to 4 inches. I tried to do this in the first place as casually as possible, working in my own office, being interrupted every five minutes or so and not remembering what I had done before the interruption. Following the discovery of microscopic particles and the development of computers, interest expanded beyond the proving ground of gravitational systems to the statistical properties of matter. In an attempt to understand the origin of irreversibility, Enrico Fermi proposed in 1953, and published in 1955, the use of the early computer MANIAC I, also at Los Alamos National Laboratory, to solve the time evolution of the equations of motion for a many body system subject to several choices of force laws. Today, this seminal work is known as the Fermi–Pasta–Ulam–Tsingou problem. The time evolution of the energy from the original work is shown in the figure to the right. In 1957, Berni Alder and Thomas Wainwright used an IBM 704 computer to simulate perfectly elastic collisions between hard spheres. In 1960, in perhaps the first realistic simulation of matter, J. B. Gibson et al. simulated radiation damage of solid copper by using a Born–Mayer type of repulsive interaction along with a cohesive surface force. In 1964, Aneesur Rahman published simulations of liquid argon that used a Lennard Jones potential; calculations of system properties, such as the coefficient of self diffusion, compared well with experimental data. Today, the Lennard Jones potential is still one of the most frequently used intermolecular potentials. It is used for describing simple substances (a. k. a. Lennard Jonesium) for conceptual and model studies and as a building block in many force fields of real substances.

Области применения и пределы

Впервые применённый в теоретической физике, метод молекулярной динамики вскоре завоевал популярность в материаловедении, а с 1970-х годов также широко используется в биохимии и биофизике. МД часто используется для уточнения трёхмерных структур белков и других макромолекул на основе экспериментальных ограничений, полученных с помощью рентгеновской кристаллографии или спектроскопии ЯМР. В физике МД применяется для изучения динамики явлений на атомном уровне, которые невозможно наблюдать непосредственно, таких как рост тонких плёнок и ионная имплантация, а также для исследования физических свойств нанотехнологических устройств, которые ещё не созданы или не могут быть созданы. В биофизике и структурной биологии метод часто используется для изучения движений макромолекул, таких как белки и нуклеиновые кислоты, что может быть полезно для интерпретации результатов определённых биофизических экспериментов и моделирования взаимодействий с другими молекулами, например, при докинге лигандов. В принципе, МД может быть использован для *ab initio* предсказания структуры белка путём моделирования сворачивания полипептидной цепи из случайной катушки. Результаты моделирования МД могут быть проверены путём сравнения с экспериментами, измеряющими молекулярную динамику, одним из популярных методов которых является спектроскопия ЯМР. Прогнозы структур, полученные с помощью МД, могут быть протестированы в рамках общедоступных экспериментов в рамках Critical Assessment of Protein Structure Prediction (CASP), хотя исторически этот метод имел ограниченный успех в этой области. Майкл Левитт, разделивший Нобелевскую премию, в том числе за применение МД к белкам, писал в 1999 году, что участники CASP обычно не использовали этот метод из-за "фундаментального недостатка молекулярной механики, а именно того, что минимизация энергии или молекулярная динамика обычно приводит к модели, менее похожей на экспериментальную структуру". Улучшение вычислительных ресурсов, позволяющее проводить более длительные траектории МД, в сочетании с современными улучшениями качества параметров силовых полей, привело к некоторым улучшениям как в предсказании структуры, так и в уточнении гомологических моделей, однако практической пользы в этих областях пока не достигнуто; многие считают параметры силовых полей ключевой областью для дальнейшей разработки. О моделировании МД сообщалось при разработке фармакофоров и проектировании лекарственных препаратов. Например, Pinto и соавторы провели моделирование МД комплексов Bcl-xL для вычисления средних положений критических аминокислот, участвующих в связывании лигандов. Carlson и соавторы использовали моделирование молекулярной динамики для идентификации соединений, которые комплементарны рецептору, вызывая при этом минимальные нарушения конформации и гибкости активного центра. Снимки белка, сделанные через равные промежутки времени в ходе моделирования, были наложены друг на друга для выявления консервативных областей связывания (консервативных как минимум в трёх из одиннадцати кадров) для разработки фармакофора. Spyrakis и соавторы использовали рабочий процесс, включающий моделирование МД, отпечатки пальцев для лигандов и белков (FLAP) и линейный дискриминантный анализ (LDA) для выявления наилучших конформаций лиганд-белка, которые могут служить шаблонами фармакофоров на основе ретроспективного ROC-анализа полученных фармакофоров. В попытке улучшить моделирование лекарственных средств на основе структуры, учитывая необходимость моделирования большого количества соединений, Hatmal и соавторы предложили комбинацию моделирования МД и анализа межмолекулярных контактов лиганд-рецептор для выявления критических межмолекулярных контактов (взаимодействий при связывании) и отсеивания избыточных в одном комплексе лиганд-белок. Критические контакты затем могут быть преобразованы в фармакофорные модели, которые могут быть использованы для виртуального скрининга. Важным фактором являются внутримолекулярные водородные связи, которые не учитываются явно в современных силовых полях, но описываются как кулоновские взаимодействия точечных зарядов атомов. Это грубое приближение, поскольку водородные связи имеют частично квантово-механическую и химическую природу. Кроме того, электростатические взаимодействия обычно рассчитываются с использованием диэлектрической проницаемости вакуума, хотя окружающий водный раствор имеет гораздо более высокую диэлектрическую проницаемость. Таким образом, использование макроскопической диэлектрической проницаемости на коротких межмолекулярных расстояниях вызывает вопросы. Наконец, ван-дер-ваальсовы взаимодействия в МД обычно описываются потенциалами Леннарда-Джонса, основанными на теории Фрица Лондона, которая применима только в вакууме. Однако все типы ван-дер-ваальсовых сил в конечном итоге имеют электростатическую природу и, следовательно, зависят от диэлектрических свойств среды. Прямое измерение сил притяжения между различными материалами (в виде постоянной Хамакера) показывает, что "взаимодействие между углеводородами через воду составляет около 10% от взаимодействия через вакуум". Для моделирования используются данные, охватывающие наносекунды (10−9 с) и микросекунды (10−6 с). Для получения этих результатов моделирования требуются дни или годы работы центрального процессора. Параллельные алгоритмы позволяют распределить нагрузку между процессорами; примером является алгоритм пространственного или силового разложения. Во время классического моделирования МД наиболее ресурсоёмкой задачей является вычисление потенциальной энергии как функции внутренних координат частиц. В рамках этого вычисления энергии наиболее дорогостоящей является не связанная или нековалентная часть. В нотации "O большое" обычные моделирования молекулярной динамики масштабируются как если все парные электростатические и ван-дер-ваальсовы взаимодействия должны учитываться явно. Эта вычислительная стоимость может быть снижена за счёт использования электростатических методов, таких как суммирование на сетке частиц Эвальда ( ), метод частиц-частиц-частиц-сетки (P3M) или методы сферического отсечения ( ). Другим фактором, влияющим на общее время работы ЦП, необходимое для моделирования, является размер временного шага интегрирования. Это промежуток времени между вычислениями потенциальной энергии. Временной шаг должен быть достаточно малым, чтобы избежать ошибок дискретизации (т.е. меньше периода, связанного с самой быстрой колебательной частотой в системе). Типичные временные шаги для классического МД составляют около 1 фемтосекунды (10−15 с). Это значение можно увеличить, используя алгоритмы, такие как алгоритм ограничения SHAKE, который фиксирует колебания самых быстрых атомов (например, водорода). Также были разработаны методы множественных временных масштабов, которые позволяют увеличить интервалы между обновлениями медленных дальнодействующих сил. Для моделирования молекул в растворителе необходимо выбрать между явным и неявным растворителем. Явные частицы растворителя (например, модели воды TIP3P, SPC/E и SPC f) должны дорогостояще рассчитываться силовым полем, в то время как неявные растворители используют подход среднего поля. Использование явного растворителя является вычислительно затратным, требуя включения примерно в десять раз большего числа частиц в моделирование. Однако гранулярность и вязкость явного растворителя необходимы для воспроизведения определённых свойств растворяемых молекул. Это особенно важно для воспроизведения химической кинетики. Во всех видах моделирования молекулярной динамики размер моделируемой ячейки должен быть достаточно большим, чтобы избежать артефактов, связанных с граничными условиями. Граничные условия часто обрабатываются путём задания фиксированных значений на границах (что может вызывать артефакты) или путём использования периодических граничных условий, при которых одна сторона моделирования замыкается на противоположную сторону, имитируя объёмную фазу (что также может вызывать артефакты).

First used in theoretical physics, the molecular dynamics method gained popularity in materials science soon afterward, and since the 1970s it has also been commonly used in biochemistry and biophysics. MD is frequently used to refine 3 dimensional structures of proteins and other macromolecules based on experimental constraints from X ray crystallography or NMR spectroscopy. In physics, MD is used to examine the dynamics of atomic level phenomena that cannot be observed directly, such as thin film growth and ion subplantation, and to examine the physical properties of nanotechnological devices that have not or cannot yet be created. In biophysics and structural biology, the method is frequently applied to study the motions of macromolecules such as proteins and nucleic acids, which can be useful for interpreting the results of certain biophysical experiments and for modeling interactions with other molecules, as in ligand docking. In principle, MD can be used for ab initio prediction of protein structure by simulating folding of the polypeptide chain from a random coil. The results of MD simulations can be tested through comparison to experiments that measure molecular dynamics, of which a popular method is NMR spectroscopy. MD derived structure predictions can be tested through community wide experiments in Critical Assessment of Protein Structure Prediction (CASP), although the method has historically had limited success in this area. Michael Levitt, who shared the Nobel Prize partly for the application of MD to proteins, wrote in 1999 that CASP participants usually did not use the method due to " a central embarrassment of molecular mechanics, namely that energy minimization or molecular dynamics generally leads to a model that is less like the experimental structure". Improvements in computational resources permitting more and longer MD trajectories, combined with modern improvements in the quality of force field parameters, have yielded some improvements in both structure prediction and homology model refinement, without reaching the point of practical utility in these areas; many identify force field parameters as a key area for further development. MD simulation has been reported for pharmacophore development and drug design. For example, Pinto et al. implemented MD simulations of Bcl xL complexes to calculate average positions of critical amino acids involved in ligand binding. Carlson et al. implemented molecular dynamics simulations to identify compounds that complement a receptor while causing minimal disruption to the conformation and flexibility of the active site. Snapshots of the protein at constant time intervals during the simulation were overlaid to identify conserved binding regions (conserved in at least three out of eleven frames) for pharmacophore development. Spyrakis et al. relied on a workflow of MD simulations, fingerprints for ligands and proteins (FLAP) and linear discriminant analysis (LDA) to identify the best ligand protein conformations to act as pharmacophore templates based on retrospective ROC analysis of the resulting pharmacophores. In an attempt to ameliorate structure based drug discovery modeling, vis à vis the need for many modeled compounds, Hatmal et al. proposed a combination of MD simulation and ligand receptor intermolecular contacts analysis to discern critical intermolecular contacts (binding interactions) from redundant ones in a single ligand–protein complex. Critical contacts can then be converted into pharmacophore models that can be used for virtual screening. An important factor is intramolecular hydrogen bonds, which are not explicitly included in modern force fields, but described as Coulomb interactions of atomic point charges. This is a crude approximation because hydrogen bonds have a partially quantum mechanical and chemical nature. Furthermore, electrostatic interactions are usually calculated using the dielectric constant of a vacuum, even though the surrounding aqueous solution has a much higher dielectric constant. Thus, using the macroscopic dielectric constant at short interatomic distances is questionable. Finally, van der Waals interactions in MD are usually described by Lennard Jones potentials based on the Fritz London theory that is only applicable in a vacuum. However, all types of van der Waals forces are ultimately of electrostatic origin and therefore depend on dielectric properties of the environment. The direct measurement of attraction forces between different materials (as Hamaker constant) shows that "the interaction between hydrocarbons across water is about 10% of that across vacuum". use data from simulations spanning nanoseconds (10−9 s) to microseconds (10−6 s). To obtain these simulations, several CPU days to CPU years are needed. Parallel algorithms allow the load to be distributed among CPUs; an example is the spatial or force decomposition algorithm. During a classical MD simulation, the most CPU intensive task is the evaluation of the potential as a function of the particles' internal coordinates. Within that energy evaluation, the most expensive one is the non bonded or non covalent part. In big O notation, common molecular dynamics simulations scale by if all pair wise electrostatic and van der Waals interactions must be accounted for explicitly. This computational cost can be reduced by employing electrostatics methods such as particle mesh Ewald summation ( ), particle particle particle mesh (P3M), or good spherical cutoff methods ( ). Another factor that impacts total CPU time needed by a simulation is the size of the integration timestep. This is the time length between evaluations of the potential. The timestep must be chosen small enough to avoid discretization errors (i. e., smaller than the period related to fastest vibrational frequency in the system). Typical timesteps for classical MD are on the order of 1 femtosecond (10−15 s). This value may be extended by using algorithms such as the SHAKE constraint algorithm, which fix the vibrations of the fastest atoms (e. g., hydrogens) into place. Multiple time scale methods have also been developed, which allow extended times between updates of slower long range forces. For simulating molecules in a solvent, a choice should be made between an explicit and implicit solvent. Explicit solvent particles (such as the TIP3P, SPC/E and SPC f water models) must be calculated expensively by the force field, while implicit solvents use a mean field approach. Using an explicit solvent is computationally expensive, requiring inclusion of roughly ten times more particles in the simulation. But the granularity and viscosity of explicit solvent is essential to reproduce certain properties of the solute molecules. This is especially important to reproduce chemical kinetics. In all kinds of molecular dynamics simulations, the simulation box size must be large enough to avoid boundary condition artifacts. Boundary conditions are often treated by choosing fixed values at the edges (which may cause artifacts), or by employing periodic boundary conditions in which one side of the simulation loops back to the opposite side, mimicking a bulk phase (which may cause artifacts too).

Микроканонический ансамбль (МКА)

В микроканоническом ансамбле система изолирована от изменений числа частиц (N), объема (V) и энергии (E). Это соответствует адиабатическому процессу без теплообмена. Траектория молекулярной динамики в микроканоническом ансамбле может рассматриваться как обмен потенциальной и кинетической энергией, при этом полная энергия сохраняется. Для системы из N частиц с координатами и скоростями , следующая пара дифференциальных уравнений первого порядка может быть записана в нотации Ньютона как

Функция потенциальной энергии системы является функцией координат частиц. В физике её часто называют просто потенциалом, а в химии – силовым полем. Первое уравнение вытекает из законов движения Ньютона; сила, действующая на каждую частицу в системе, может быть рассчитана как отрицательный градиент . Для каждого шага по времени положение и скорость каждой частицы могут быть проинтегрированы с использованием метода симплектического интегрирования, такого как интегрирование Верле. Временная эволюция и называется траекторией. Зная начальные положения (например, полученные из теоретических соображений) и скорости (например, случайные, распределенные по Гауссу), можно рассчитать все будущие (или прошлые) положения и скорости. Одним из частых источников путаницы является значение температуры в молекулярной динамике. Обычно мы имеем дело с макроскопическими температурами, которые характеризуют огромное количество частиц, но температура – это статистическая величина. Если число атомов достаточно велико, статистическую температуру можно оценить на основе мгновенной температуры, которая определяется путем приравнивания кинетической энергии системы к nkBT/2, где n – число степеней свободы системы. Явление, связанное с температурой, возникает из-за небольшого числа атомов, используемых в расчетах молекулярной динамики. Например, рассмотрим моделирование роста медной пленки, начиная с подложки, содержащей 500 атомов, и энергией осаждения 100 эВ. В реальном мире 100 эВ от осажденного атома быстро передавались бы и распределялись между большим числом атомов (или более) без существенного изменения температуры. Однако, когда в системе всего 500 атомов, подложка почти мгновенно испаряется под воздействием осаждения. Подобное происходит и в биофизических симуляциях. Температура системы в NVE естественным образом повышается, когда макромолекулы, такие как белки, претерпевают экзотермические конформационные изменения и образуют связи.

Канонический ансамбль (NVT)

В каноническом ансамбле сохраняются количество вещества (N), объем (V) и температура (T). Его также иногда называют молекулярной динамикой с постоянной температурой (CTMD). В ансамбле NVT энергия эндотермических и экзотермических процессов обменивается с термостатом. Существует множество алгоритмов термостатирования, позволяющих добавлять и удалять энергию из границ MD-симуляции более или менее реалистичным образом, приближая канонический ансамбль. Популярные методы контроля температуры включают в себя масштабирование скоростей, термостат Нозе-Гувера, цепи Нозе-Гувера, термостат Берендсена, термостат Андерсена и динамику Ланжевена. Термостат Берендсена может приводить к эффекту "летающего кубика льда", вызывающему нефизические трансляции и вращения моделируемой системы. Получение канонического распределения конформаций и скоростей с использованием этих алгоритмов – нетривиальная задача. Зависимость этого результата от размера системы, выбора термостата, его параметров, шага по времени и используемого интегратора является предметом многочисленных исследований в данной области.

Изотермальный изобарический (NPT) ансамбль

В изотермическо-изобарическом ансамбле сохраняются количество вещества (N), давление (P) и температура (T). Помимо термостата, необходим баростат. Это наиболее близко соответствует лабораторным условиям с открытой колбой, находящейся в контакте с температурой и давлением окружающей среды. При моделировании биологических мембран изотропный контроль давления неприемлем. Для липидных бислоев контроль давления осуществляется при постоянной площади мембраны (NPAT) или постоянном поверхностном натяжении "гамма" (NPγT).

Обобщенные ансамбли

Метод обмена репликами является обобщенным ансамблем. Он был первоначально разработан для преодоления медленной динамики в неупорядоченных спиновых системах. Он также известен как параллельный отжиг. В формулировке молекулярной динамики с обменом репликами (REMD) предпринимается попытка решить проблему множественных минимумов путем обмена температурой между не взаимодействующими репликами системы, работающими при различных температурах.

Потенциалы в симуляциях МД

Молекулярное моделирование требует определения потенциальной функции или описания членов, посредством которых частицы в моделировании будут взаимодействовать. В химии и биологии это обычно называют силовым полем, а в физике материалов – межатомным потенциалом. Потенциалы могут быть определены на различных уровнях физической точности; наиболее часто используемые в химии основаны на молекулярной механике и включают в себя классическое описание взаимодействия частиц, способное воспроизводить структурные и конформационные изменения, но обычно не воспроизводящее химические реакции. Переход от полностью квантового описания к классическому потенциалу предполагает два основных приближения. Первое – это приближение Борна-Оппенгеймера, которое утверждает, что динамика электронов настолько быстра, что они могут считаться мгновенно реагирующими на движение ядер. Как следствие, их можно рассматривать отдельно. Второе приближение рассматривает ядра, которые значительно тяжелее электронов, как точечные частицы, подчиняющиеся классической ньютоновской динамике. В классической молекулярной динамике влияние электронов аппроксимируется одной поверхностью потенциальной энергии, обычно представляющей основное состояние. Когда требуется более высокая детализация, используются потенциалы, основанные на квантовой механике; некоторые методы стремятся создать гибридные классические/квантовые потенциалы, в которых основная часть системы описывается классически, а небольшая область – как квантовая система, обычно претерпевающая химическую трансформацию.

Эмпирические потенциалы

Эмпирические потенциалы, используемые в химии, часто называются силовыми полями, а те, что используются в физике материалов, – межатомными потенциалами. Большинство силовых полей в химии являются эмпирическими и состоят из суммы связанных сил, соответствующих химическим связям, углам между связями и торсионным углам, и несвязанных сил, соответствующих силам Ван-дер-Ваальса и электростатическому заряду. Эмпирические потенциалы представляют квантово-механические эффекты в ограниченном виде посредством ad hoc функциональных приближений. Эти потенциалы содержат свободные параметры, такие как атомный заряд, параметры Ван-дер-Ваальса, отражающие оценки атомного радиуса, и равновесную длину связи, угол и торсионный угол; они определяются путем подгонки к детальным электронным расчетам (квантово-химическим симуляциям) или экспериментальным физическим свойствам, таким как упругие константы, параметры решетки и спектроскопические измерения. Из-за нелокального характера несвязанных взаимодействий, они включают как минимум слабые взаимодействия между всеми частицами в системе. Их вычисление обычно является узким местом, ограничивающим скорость MD-моделирования. Для снижения вычислительных затрат силовые поля используют численные приближения, такие как смещенные радиусы отсечения, алгоритмы реакционного поля, суммирование ячечных частиц Эвальда или более современный метод частиц-частиц-частиц-ячеек (P3M). Химические силовые поля обычно используют заданные схемы связывания (за исключением динамики ab initio) и, следовательно, не могут явно моделировать процесс разрыва химических связей и протекание реакций. С другой стороны, многие потенциалы, используемые в физике, такие как основанные на формализме порядка связей, могут описывать различные координации системы и разрыв связей. Примерами таких потенциалов являются потенциал Бреннера для углеводородов и его дальнейшие разработки для систем C Si H и C O H. Потенциал ReaxFF можно рассматривать как полностью реактивный гибрид между потенциалами порядка связей и химическими силовыми полями.

Потенциалы пар против потенциалов многих тел

Потенциальные функции, представляющие несвязанную энергию, формулируются как сумма взаимодействий между частицами системы. Самый простой выбор, используемый во многих популярных силовых полях, – это "парный потенциал", в котором общая потенциальная энергия может быть рассчитана как сумма энергетических вкладов между парами атомов. Поэтому эти силовые поля также называются "аддитивными силовыми полями". Примером такого парного потенциала является не связанный потенциал Леннарда-Джонса (также называемый потенциалом 6–12), используемый для расчета сил Ван-дер-Ваальса. Другим примером является модель ионной решетки Борна (ионная). Первый член в следующем уравнении – закон Кулона для пары ионов, второй член – отталкивание на коротких расстояниях, объясняемое принципом исключения Паули, а последний член – член дисперсионного взаимодействия. Обычно моделирование включает только дипольный член, хотя иногда также включается квадрапольный член. При nl = 6 этот потенциал также называется потенциалом Кулона–Бакингема. В многочастичных потенциалах потенциальная энергия включает в себя влияние трех и более взаимодействующих друг с другом частиц. В симуляциях с парными потенциалами глобальные взаимодействия в системе также существуют, но они возникают только через парные члены. В многочастичных потенциалах потенциальную энергию нельзя определить как сумму по парам атомов, поскольку эти взаимодействия рассчитываются явно как комбинация членов более высокого порядка. С точки зрения статистики, зависимость между переменными в общем случае не может быть выражена только парными произведениями степеней свободы. Например, потенциал Терсоффа, который первоначально использовался для моделирования углерода, кремния и германия, а затем и для широкого спектра других материалов, включает в себя сумму по группам из трех атомов, при этом углы между атомами являются важным фактором потенциала. Другими примерами являются метод встроенных атомов (EAM) и EDIP, где электронная плотность состояний в области атома рассчитывается как сумма вкладов от окружающих атомов, а вклад потенциальной энергии является функцией этой суммы.

The potential functions representing the non bonded energy are formulated as a sum over interactions between the particles of the system. The simplest choice, employed in many popular force fields, is the "pair potential", in which the total potential energy can be calculated from the sum of energy contributions between pairs of atoms. Therefore, these force fields are also called "additive force fields". An example of such a pair potential is the non bonded Lennard Jones potential (also termed the 6–12 potential), used for calculating van der Waals forces. Another example is the Born (ionic) model of the ionic lattice. The first term in the next equation is Coulomb's law for a pair of ions, the second term is the short range repulsion explained by Pauli's exclusion principle and the final term is the dispersion interaction term. Usually, a simulation only includes the dipolar term, although sometimes the quadrupolar term is also included. When nl = 6, this potential is also called the Coulomb–Buckingham potential. In many body potentials, the potential energy includes the effects of three or more particles interacting with each other. In simulations with pairwise potentials, global interactions in the system also exist, but they occur only through pairwise terms. In many body potentials, the potential energy cannot be found by a sum over pairs of atoms, as these interactions are calculated explicitly as a combination of higher order terms. In the statistical view, the dependency between the variables cannot in general be expressed using only pairwise products of the degrees of freedom. For example, the Tersoff potential, which was originally used to simulate carbon, silicon, and germanium, and has since been used for a wide range of other materials, involves a sum over groups of three atoms, with the angles between the atoms being an important factor in the potential. Other examples are the embedded atom method (EAM), the EDIP, where the electron density of states in the region of an atom is calculated from a sum of contributions from surrounding atoms, and the potential energy contribution is then a function of this sum.

Полуэмпирические потенциалы

Полуэмпирические потенциалы используют матричное представление из квантовой механики. Однако значения матричных элементов определяются с помощью эмпирических формул, оценивающих степень перекрытия конкретных атомных орбиталей. Затем матрица диагонализируется для определения заполнения различных атомных орбиталей, и эмпирические формулы вновь используются для определения энергетического вклада этих орбиталей. Существует большое разнообразие полуэмпирических потенциалов, называемых потенциалами жесткого связывания, которые различаются в зависимости от моделируемых атомов.

Поляризуемые потенциалы

Большинство классических силовых полей неявно учитывают эффект поляризуемости, например, путем масштабирования частичных зарядов, полученных из квантово-химических расчетов. Эти частичные заряды фиксированы относительно массы атома. Однако, моделирование молекулярной динамики может явно моделировать поляризуемость, вводя индуцированные диполи различными методами, такими как частицы Друде или флуктуирующие заряды. Это позволяет динамически перераспределять заряд между атомами в ответ на локальное химическое окружение. На протяжении многих лет поляризуемые MD-симуляции считались следующим поколением методов. Для однородных жидкостей, таких как вода, повышение точности было достигнуто благодаря включению поляризуемости. Также получены многообещающие результаты для белков. Тем не менее, до сих пор неясно, как наилучшим образом аппроксимировать поляризуемость в симуляции. Этот вопрос становится особенно важным, когда частица испытывает различные среды в ходе своей траектории моделирования, например, при транслокации лекарственного препарата через клеточную мембрану.

Потенциалы в методах ab initio

В классической молекулярной динамике одна потенциальная энергетическая поверхность (обычно основное состояние) описывается силовым полем. Это является следствием приближения Борна — Оппенгеймера. В возбужденных состояниях, при химических реакциях или когда требуется более точное описание, электронное поведение можно получить из первых принципов, используя квантово-механический метод, такой как теория функционала плотности. Это называется Ab Initio Molecular Dynamics (AIMD). Из-за вычислительных затрат на обработку электронных степеней свободы, вычислительная нагрузка при таких симуляциях значительно выше, чем в классической молекулярной динамике. По этой причине AIMD обычно ограничивается меньшими системами и более короткими временными интервалами. Для расчета потенциальной энергии системы «на лету», необходимого для определения конформаций в траектории, могут использоваться квантово-механические и химические методы. Этот расчет обычно выполняется вблизи реакционной координаты. Хотя могут применяться различные приближения, они основаны на теоретических соображениях, а не на эмпирической подгонке. Расчеты ab initio предоставляют большой объем информации, недоступной при использовании эмпирических методов, например, плотность электронных состояний или другие электронные свойства. Важным преимуществом использования методов ab initio является возможность изучения реакций, связанных с разрывом или образованием ковалентных связей, которые соответствуют нескольким электронным состояниям. Более того, методы ab initio позволяют учитывать эффекты, выходящие за рамки приближения Борна — Оппенгеймера, используя подходы, такие как смешанная квантово-классическая динамика.

Гибридный QM/MM

Методы КМ (квантовой механики) очень мощные. Однако они вычислительно затратны, в то время как методы MM (классической или молекулярной механики) быстры, но имеют ряд ограничений (требуют обширной параметризации; полученные оценки энергии недостаточно точны; не могут использоваться для моделирования реакций, в которых разрываются/образуются ковалентные связи; и ограничены в своей способности предоставлять точные детали об химическом окружении). Появился новый класс методов, сочетающий преимущества расчетов КМ (точность) и MM (скорость). Эти методы называются смешанными или гибридными квантово-механическими и молекулярно-механическими методами (гибридные QM/MM). Наиболее важное преимущество гибридного метода QM/MM – это скорость. Вычислительная стоимость классической молекулярной динамики (МД) в простейшем случае масштабируется как O(n²), где n – число атомов в системе. Это в основном связано с членом, описывающим электростатические взаимодействия (каждая частица взаимодействует с каждой другой частицей). Однако использование радиуса отсечения, периодического обновления списков пар и, в последнее время, модификаций метода частиц в ячейке Эвальда (PME) снизило эту зависимость до O(n) – O(n²). Иными словами, моделирование системы с вдвое большим числом атомов потребует от двух до четырех раз больше вычислительных ресурсов. С другой стороны, простейшие ab initio расчеты обычно масштабируются как O(n³) или хуже (для ограниченных расчетов Хартри-Фока предложена зависимость ~O(n²․⁷)). Чтобы преодолеть это ограничение, небольшая часть системы обрабатывается квантово-механически (обычно активный центр фермента), а остальная часть – классически. В более сложных реализациях методы QM/MM позволяют учитывать как легкие ядра, подверженные квантовым эффектам (например, водород), так и электронные состояния. Это позволяет генерировать волновые функции водорода (аналогичные электронным волновым функциям). Эта методология оказалась полезной при изучении таких явлений, как туннелирование водорода. Одним из примеров, где методы QM/MM привели к новым открытиям, является расчет переноса гидрида в ферменте печеночной алкогольдегидрогеназе. В этом случае квантовое туннелирование играет важную роль для водорода, поскольку оно определяет скорость реакции.

QM (quantum mechanical) methods are very powerful. However, they are computationally expensive, while the MM (classical or molecular mechanics) methods are fast but suffer from several limits (require extensive parameterization; energy estimates obtained are not very accurate; cannot be used to simulate reactions where covalent bonds are broken/formed; and are limited in their abilities for providing accurate details regarding the chemical environment). A new class of method has emerged that combines the good points of QM (accuracy) and MM (speed) calculations. These methods are termed mixed or hybrid quantum mechanical and molecular mechanics methods (hybrid QM/MM). The most important advantage of hybrid QM/MM method is the speed. The cost of doing classical molecular dynamics (MM) in the most straightforward case scales O(n2), where n is the number of atoms in the system. This is mainly due to electrostatic interactions term (every particle interacts with every other particle). However, use of cutoff radius, periodic pair list updates and more recently the variations of the particle mesh Ewald's (PME) method has reduced this to between O(n) to O(n2). In other words, if a system with twice as many atoms is simulated then it would take between two and four times as much computing power. On the other hand, the simplest ab initio calculations typically scale O(n3) or worse (restricted Hartree–Fock calculations have been suggested to scale ~O(n2.7)). To overcome the limit, a small part of the system is treated quantum mechanically (typically active site of an enzyme) and the remaining system is treated classically. In more sophisticated implementations, QM/MM methods exist to treat both light nuclei susceptible to quantum effects (such as hydrogens) and electronic states. This allows generating hydrogen wave functions (similar to electronic wave functions). This methodology has been useful in investigating phenomena such as hydrogen tunneling. One example where QM/MM methods have provided new discoveries is the calculation of hydride transfer in the enzyme liver alcohol dehydrogenase. In this case, quantum tunneling is important for the hydrogen, as it determines the reaction rate.

Силы дальнего действия

Взаимодействие на большом расстоянии – это взаимодействие, в котором пространственное затухание взаимодействия не быстрее, чем , где – размерность системы. Примеры включают электростатическое взаимодействие между ионами и диполь-дипольное взаимодействие между молекулами. Моделирование этих сил представляет собой значительную сложность, поскольку они существенны на расстояниях, которые могут превышать половину длины расчетной области при моделировании систем, содержащих тысячи частиц. Хотя одним из решений было бы значительное увеличение длины расчетной области, этот прямой подход не является оптимальным, так как вычислительные затраты на моделирование существенно возрастут. Сферическое усечение потенциала также неприемлемо, поскольку при приближении к расстоянию отсечения могут наблюдаться нереалистичные эффекты.

Управляемая молекулярная динамика (SMD)

Симуляции направленной молекулярной динамики (SMD), или симуляции силового зонда, применяют силы к белку для манипулирования его структурой, растягивая его вдоль желаемых степеней свободы. Эти эксперименты позволяют выявить структурные изменения в белке на атомном уровне. SMD часто используется для моделирования событий, таких как механическое разворачивание или растяжение. Существуют два типичных протокола SMD: один, в котором поддерживается постоянная скорость растяжения, и один, в котором поддерживается постоянная приложенная сила. Обычно часть исследуемой системы (например, атом в белке) фиксируется гармоническим потенциалом. Затем к определенным атомам прикладываются силы с постоянной скоростью или постоянной силой. Метод зонтичной выборки используется для перемещения системы вдоль желаемой реакционной координаты путем изменения, например, сил, расстояний и углов, которыми манипулируют в симуляции. Благодаря зонтичной выборке адекватно исследуются все конфигурации системы – как высокоэнергетические, так и низкоэнергетические. Затем изменение свободной энергии каждой конфигурации может быть рассчитано как потенциал средней силы (PMF). Популярным методом вычисления PMF является метод анализа взвешенных гистограмм (WHAM), который анализирует серию симуляций зонтичной выборки. Многие важные применения SMD находятся в области разработки лекарств и биомолекулярных наук. Например, SMD использовался для исследования стабильности протофибрилл болезни Альцгеймера, для изучения взаимодействия белок-лиганд в циклинозависимой киназе 5 и даже для демонстрации влияния электрического поля на комплекс тромбина (белка) и аптамера (нуклеотида) среди множества других интересных исследований.

Примеры применения

Молекулярная динамика используется во многих областях науки. Первое моделирование упрощенного биологического процесса сворачивания было опубликовано в 1975 году. Его моделирование, опубликованное в Nature, проложило путь к обширной области современного вычислительного сворачивания белков. Первая симуляция биологического процесса была опубликована в 1976 году. Его моделирование, опубликованное в Nature, проложило путь к пониманию движения белка как существенного для функции, а не просто вспомогательного. MD – стандартный метод для обработки каскадов столкновений в режиме теплового всплеска, то есть эффектов, которые энергетическое нейтронное и ионное облучение оказывает на твердые тела и поверхности. Следующие биофизические примеры иллюстрируют значительные усилия по созданию симуляций систем очень больших размеров (полный вирус) или очень длительного времени симуляции (до 1,112 миллисекунд):
Симуляция методом MD полного спутникового вируса табачной мозаики (STMV) (2006, Размер: 1 миллион атомов, Время симуляции: 50 нс, Программа: NAMD). Этот вирус представляет собой небольшой икосаэдрический вирус растений, который усугубляет симптомы инфекции вирусом табачной мозаики (TMV). Молекулярная динамика была использована для изучения механизмов сборки вируса. Вся частица STMV состоит из 60 идентичных копий одного белка, составляющего вирусную капсиду (оболочку), и 1063-нуклеотидного одноцепочечного РНК-генома. Ключевой вывод заключается в том, что капсид очень нестабилен в отсутствие РНК внутри. Для выполнения симуляции одному настольному компьютеру 2006 года потребовалось бы около 35 лет. Поэтому она была выполнена на множестве процессоров параллельно с непрерывным обменом данными между ними. Симуляции сворачивания головки виллина в деталях на уровне атома (2006, Размер: 20 000 атомов; Время симуляции: 500 мкс = 500 000 нс, Программа: Folding@home). Эта симуляция была выполнена на 200 000 процессорах персональных компьютеров по всему миру. На этих компьютерах была установлена программа Folding@home – крупномасштабный распределенный вычислительный проект, координируемый Виджаем Панде в Стэнфордском университете. Кинетические свойства белка Villin Headpiece были исследованы с использованием множества независимых коротких траекторий, запущенных процессорами без непрерывной связи в реальном времени. Одним из использованных методов был анализ значения Pfold, который измеряет вероятность сворачивания перед разворачиванием конкретной исходной конформации. Pfold предоставляет информацию о структурах переходного состояния и упорядочении конформаций вдоль пути сворачивания. Каждая траектория в расчете Pfold может быть относительно короткой, но требуется множество независимых траекторий. Длительные непрерывные траектории симуляции были выполнены на Anton – массивно параллельном суперкомпьютере, разработанном и построенном на основе специализированных интегральных схем (ASIC) и межсоединений компанией D. E. Shaw Research. Самый длинный опубликованный результат симуляции, выполненной с использованием Anton, – это 1,112-миллисекундная симуляция NTL9 при 355 K; также была выполнена вторая, независимая 1,073-миллисекундная симуляция этой конфигурации (и множество других симуляций с непрерывным химическим временем более 250 мкс). В книге «Как быстро сворачиваются белки» исследователи Крестен Линдорф Ларсен, Стефано Пиана, Рон О. Дрор и Дэвид Э. Шоу обсуждают «результаты моделирования молекулярной динамики на атомном уровне в течение периодов от 100 мкс до 1 мс, которые выявляют набор общих принципов, лежащих в основе сворачивания 12 структурно различных белков». Изучение этих разнообразных длинных траекторий, обеспеченное специализированным пользовательским оборудованием, позволяет им заключить, что «в большинстве случаев сворачивание следует по одному доминирующему маршруту, в котором элементы нативной структуры появляются в порядке, сильно коррелирующем с их склонностью к формированию в развернутом состоянии». Другое важное применение метода MD выигрывает от его способности к трехмерной характеристике и анализу микроструктурной эволюции в атомном масштабе. Симуляции MD используются для характеристики эволюции размера зерна, например, при описании износа и трения нанокристаллических материалов Al и Al(Zr). Эволюция дислокаций и эволюция размера зерна анализируются во время процесса трения в этой симуляции. Поскольку метод MD предоставил полную информацию о микроструктуре, эволюция размера зерна была рассчитана в 3D с использованием методов сопоставления полиэдральных шаблонов, сегментации зерен и кластеризации графов. В такой симуляции метод MD обеспечил точное измерение размера зерна. Используя эту информацию, были извлечены, измерены и представлены фактические структуры зерен. По сравнению с традиционным методом использования СЭМ с одним двумерным срезом материала, MD обеспечивает трехмерный и точный способ характеристики микроструктурной эволюции в атомном масштабе.

Molecular dynamics is used in many fields of science. First MD simulation of a simplified biological folding process was published in 1975. Its simulation published in Nature paved the way for the vast area of modern computational protein folding. First MD simulation of a biological process was published in 1976. Its simulation published in Nature paved the way for understanding protein motion as essential in function and not just accessory. MD is the standard method to treat collision cascades in the heat spike regime, i. e., the effects that energetic neutron and ion irradiation have on solids and solid surfaces. The following biophysical examples illustrate notable efforts to produce simulations of a systems of very large size (a complete virus) or very long simulation times (up to 1.112 milliseconds):
MD simulation of the full satellite tobacco mosaic virus (STMV) (2006, Size: 1 million atoms, Simulation time: 50 ns, program: NAMD) This virus is a small, icosahedral plant virus that worsens the symptoms of infection by Tobacco Mosaic Virus (TMV). Molecular dynamics simulations were used to probe the mechanisms of viral assembly. The entire STMV particle consists of 60 identical copies of one protein that make up the viral capsid (coating), and a 1063 nucleotide single stranded RNA genome. One key finding is that the capsid is very unstable when there is no RNA inside. The simulation would take one 2006 desktop computer around 35 years to complete. It was thus done in many processors in parallel with continuous communication between them. Folding simulations of the Villin Headpiece in all atom detail (2006, Size: 20,000 atoms; Simulation time: 500 μs= 500,000 ns, Program: Folding@home) This simulation was run in 200,000 CPU's of participating personal computers around the world. These computers had the Folding@home program installed, a large scale distributed computing effort coordinated by Vijay Pande at Stanford University. The kinetic properties of the Villin Headpiece protein were probed by using many independent, short trajectories run by CPU's without continuous real time communication. One method employed was the Pfold value analysis, which measures the probability of folding before unfolding of a specific starting conformation. Pfold gives information about transition state structures and an ordering of conformations along the folding pathway. Each trajectory in a Pfold calculation can be relatively short, but many independent trajectories are needed. Long continuous trajectory simulations have been performed on Anton, a massively parallel supercomputer designed and built around custom application specific integrated circuits (ASICs) and interconnects by D. E. Shaw Research. The longest published result of a simulation performed using Anton is a 1.112 millisecond simulation of NTL9 at 355 K; a second, independent 1.073 millisecond simulation of this configuration was also performed (and many other simulations of over 250 μs continuous chemical time). In How Fast Folding Proteins Fold, researchers Kresten Lindorff Larsen, Stefano Piana, Ron O. Dror, and David E. Shaw discuss "the results of atomic level molecular dynamics simulations, over periods ranging between 100 μs and 1 ms, that reveal a set of common principles underlying the folding of 12 structurally diverse proteins." Examination of these diverse long trajectories, enabled by specialized, custom hardware, allow them to conclude that "In most cases, folding follows a single dominant route in which elements of the native structure appear in an order highly correlated with their propensity to form in the unfolded state." Another important application of MD method benefits from its ability of 3 dimensional characterization and analysis of microstructural evolution at atomic scale. MD simulations are used in characterization of grain size evolution, for example, when describing wear and friction of nanocrystalline Al and Al(Zr) materials. Dislocations evolution and grain size evolution are analyzed during the friction process in this simulation. Since MD method provided the full information of the microstructure, the grain size evolution was calculated in 3D using the Polyhedral Template Matching, Grain Segmentation, and Graph clustering methods. In such simulation, MD method provided an accurate measurement of grain size. Making use of these information, the actual grain structures were extracted, measured, and presented. Compared to the traditional method of using SEM with a single 2 dimensional slice of the material, MD provides a 3 dimensional and accurate way to characterize the microstructural evolution at atomic scale.