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

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

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

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

Квантовая механика и уравнение Шрёдингера

Мы начнем с квантовой механики в упрощенной ситуации, когда изучаемая нами система находится в равновесии. Это описывается не зависящим от времени уравнением Шрёдингера:

Здесь m — масса, p — импульс, V — потенциальная энергия, E — полная энергия системы.

Мы хотим решить это уравнение для ψ, которая называется волновой функцией. Он кодирует плотность изучаемых нами частиц в каждой точке пространства.

Хотя на первый взгляд это может показаться довольно коротким и простым, за этим скрывается много сложностей, как концептуальных, так и вычислительных. В этом блоге я сосредоточусь на вычислительной сложности.

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

Импульс p на самом деле является производной по положению в квантовой механике. Точнее, p = — i ℏ ∂, где ℏ — очень малое число, называемое постоянной Планка. Его малость — причина, по которой мы можем проживать повседневную жизнь, не думая о квантовой механике.

Это приводит к чему-то более проблематичному, что я также заместил под ковер. Поскольку импульс является производной по положению, действующему на волновую функцию, а каждая частица имеет свой собственный импульс, волновая функция является не просто функцией положения, а функцией положений каждого одна частица. Блок материала размером с ваш смартфон имеет около 10²⁵ электронов и ядер, это много аргументов!

Между прочим, потенциал — это кулоновский потенциал, с которым заряженные частицы действуют друг на друга, так что он также зависит от всех координат.

Чтобы лучше понять вычислительную сложность, полезно посмотреть, как мы будем решать это численно. По сути, большинство численных методов решают дифференциальное уравнение, оценивая само уравнение по определенному количеству точек. Это позволяет преобразовать дифференциальное уравнение в линейную алгебру, матричное уравнение по существу вида A x = b, где A — матрица, а b — вектор, которые мы оба получаем из этой процедуры, и если мы решим это для вектора x, мы получим аппроксимация решения, которое мы ищем, на том же наборе точек.

Сколько точек нам нужно, сильно зависит от деталей, но давайте для простоты скажем, что система ограничена одной линией, другими словами, пространство только одномерное. Давайте также ограничим его интервалом между 0 и 1 и скажем, что мы используем K точек, расположенных в позициях x_i = i /K на интервале.

Возможно, K=20 может быть достаточно в какой-то очень простой ситуации, так что это не так уж плохо. Помните, однако, что каждая отдельная частица имеет положение, и волновая функция зависит от всех них. Итак, нам нужны все комбинации позиций. Таким образом, для системы с 10²⁵ частиц это составляет 20 баллов в степени 10²⁵.

Другой способ думать об этом состоит в том, что если мы добавим в нашу систему только один электрон, нам теперь нужно будет решить исходную систему один раз для каждого возможного местоположения добавленного электрона. Так что в нашем примере это сделало бы проблему в 20 раз больше, это известно как экспоненциальная сложность.

Эти нелепые цифры получены из сочетания экспоненциальной сложности проблемы с уже огромным количеством задействованных частиц.

Это означает, что прямое использование этого уравнения возможно только для самых маленьких систем, что на самом деле означает практически только атом водорода: одно ядро ​​с одним электроном.

приближения

Ясно, что тогда нам нужно сделать некоторые приближения к этому уравнению. Я концептуально проиллюстрирую серию очень стандартных объяснений.

Первый и простейший шаг следует из наблюдения, что протон намного тяжелее электрона примерно в 1000 раз. Так что даже для водорода, ядро ​​которого состоит из одного протона, это существенный фактор, а для ядер большего размера соотношение еще больше, например, углерод содержит 12 протонов и нейтронов (которые имеют аналогичную массу протонов), поэтому их масса примерно в 10⁴ раз больше, чем у электрона.

Теперь тяжелые объекты имеют тенденцию двигаться медленнее, чем легкие, в этом суть того, что называется приближение зажатых ядер.

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

И наоборот, с точки зрения ядра, к тому времени, когда оно само движется, у всех электронов будет достаточно времени, чтобы успокоиться в своих равновесных положениях, в так называемом их основном состоянии.

Так что это позволяет нам в некотором смысле отделить электроны от ядер. Само по себе это пока не дает значительного улучшения, поэтому нужно упрощать каждую часть по отдельности.

Электроны: теория функционала плотности (DFT)

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

Что мы можем сделать, так это упростить все взаимодействия до взаимодействий между каждой отдельной частицей и некоторой усредненной частицей, это известно как приближение среднего поля. То есть мы представляем, что все электроны вносят вклад в среднюю электронную плотность. Эта плотность создает некоторую потенциальную энергию, которую ощущает каждый электрон. Если это звучит рекурсивно, то это потому, что так оно и есть. Это цена, которую мы платим за такое приближение: мы получаем более простую задачу, но нам приходится решать ее много раз, пока она не станет самосогласованной.

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

Этот подход называется теорией функционала плотности (DFT).

Что я упустил, так это то, что представляет собой этот потенциал, который ощущает каждый электрон из-за общей плотности. Это то, что дало название этой теории: функционал плотности. Функционал — это то, что математики называют функцией, аргумент которой сам является функцией. В этом случае потенциал является функцией электронной плотности, которая сама является функцией положения (на этот раз одного положения!).

Интересная ситуация заключается в том, что все это «приближение» известно как точное. Так что на данный момент это вовсе не приближение, мы не теряем точности при этом. Однако есть одна загвоздка. Теория только говорит нам, что соответствующий функционал плотности существует, что делает его точным в любой ситуации, но не говорит нам что это .

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

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

Конкретно, в примере, где мы использовали 10 узлов сетки, переход от 100 к 101 частице с использованием уравнения Шредингера напрямую в 10 раз медленнее, тогда как DFT становится только на 3% медленнее. Вы видите, насколько важно это масштабирование вычислительной сложности.

На практике, в то время как уравнение Шредингера можно было применить почти только к простейшим атомам, DFT можно использовать не только для отдельных атомов, но и для молекул и даже для нескольких молекул вместе. Тем не менее, это все еще очень дорого в вычислительном отношении, это не то, что вы можете запустить на своем ноутбуке, а скорее на кластере, и, возможно, вам придется ждать несколько недель в зависимости от проблемы. Хотя с того, с чего мы начали, это уже большой шаг вперед.

Ядра: молекулярная динамика (МД)

Теперь вернемся к ядрам. Здесь наиболее важным источником приближения является их масса. Мы уже использовали это, чтобы отделить их от электронов, так как их масса намного больше. Теперь мы можем использовать это по-другому.

Квантовая механика становится необходимой для описания системы только в том случае, если она очень мала и легка. Для движения бильярдных шаров по бильярдному столу мы, к счастью, можем просто использовать классическую механику. Формально это предел квантовой механики, прозаически известный как классический предел. Помните, что в уравнении Шредингера я сказал, что импульс на самом деле является производной по положению, а волновая функция кодирует вероятность найти где-то частицу. В классическом пределе импульс снова становится просто числом, а вероятность исчезает, так что частица просто описывается ее положением и импульсом.

Удивительно, но даже для чего-то такого маленького, как ядро, это хорошее приближение. Это верно для всех ядер, кроме самых маленьких, таких как водород.

Это приводит к очень простой системе, в которой у нас есть потенциал, исходящий от электронов в дополнение к кулоновскому потенциалу между самими ядрами. Мы получаем силы, просто взяв производную от потенциала, и с помощью этих сил мы можем развивать ядра во времени.

Этот подход называется Молекулярная динамика (МД) и показан на блок-схеме ниже.

Здесь снова пришлось ввести неизвестный потенциал, в данном случае названный межатомным потенциалом. Это фиксирует влияние электронов на ядра. И снова это априори неизвестно, поэтому приходится догадываться.

В этом случае есть очень естественный кандидат. Мы только что обсудили способ вычисления электронной плотности для любой конфигурации ядер, а именно DFT. Так что мы могли бы сделать это, а затем просто использовать кулоновский потенциал, исходящий из этой электронной плотности. Это известно как молекулярная динамика ab-initio (aiMD), потому что нам не нужно было делать еще одно предположение о потенциале.

Преимущество aiMD в том, что он обычно очень точен. Однако он требует полного вычисления ДПФ для каждого шага МД, что делает его очень дорогим в вычислительном отношении.

Альтернатива состоит в том, чтобы придумать простую форму межатомного потенциала, которую можно очень быстро оценить и которая может быть точной для определенных систем. Таких догадок существует множество.

Так что здесь вычислительная сложность очень сильно зависит от используемого межатомного потенциала, остальные шаги обычно незначительны по сравнению с ним. Менее дорогие потенциалы можно использовать для моделирования миллионов атомов, например, в больших белках или кристаллических структурах.

Краткое содержание

Диаграмма ниже суммирует то, что мы обсуждали.
С помощью ряда физически мотивированных и во многих ситуациях очень точных приближений мы сократили элегантную, но трудноразрешимую в вычислительном отношении формулировку уравнения Шредингера на две части: теорию функционала плотности (DFT) для электронов и молекулярную динамику (MD ) для ядер.

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

Именно здесь на сцену выходит машинное обучение: вместо того, чтобы самим угадывать эти потенциалы, мы можем позволить модели машинного обучения изучить их на основе данных. Мы обсудим это во второй части этой серии.