Расчет структуры молекул по программам молекулярная механика. Эмпирические методы молекулярной механики. Достоинства и недостатки методов молекулярной механики

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

Простейшие модели молекулярной механики учитывают растяжения связей, деформацию валентных и двугранных (торсионных) углов, взаимодействие валентно несвязанных атомов, называемое также ван-дер-ваальсовым взаимодействием, электростатические вклады и т.д.:

U = U раст + U деф + U торс + U вдв + U эл-стат.

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

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

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

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

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

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

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

Для реализации квантовомеханического подхода к исследованию строения молекул необходимо сделать ряд упрощений. С этой целью чаще всего используют приближение самосогласованного поля Хартри-Фока. Точность расчетов «из первых принципов» для геометрического строения молекул изменяется в достаточно широких пределах (от умеренной до очень высокой) в зависимости от вида применяемых атомных волновых функций. Так, если вычисления проводятся с расширенным базисным набором, включающим ((-орбитали для элементов второго периода и р-орбитали атома водорода, то структурные и энергетические данные получаются с точностью, сравнимой с результатами молекулярной механики, в которой предел точности обусловлен не методом расчета, а разницей в строгом физическом определении свойств молекул, например разницей в длинах связей. Следует отметить, что подробные квантовомеханические расчеты геометрического строения молекул можно выполнить для любой молекулы или любого фрагмента без привлечения каких-либо экспериментальных сведений об изучаемой системе.

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

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

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

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

Если молекула, состоящая из N атомов и описываемая ЗгУ координатами х 1У деформируется по отношению к своей равновесной конфигурации с энергией?/ 0 и координатами х 0 , то ее потенциальная энергия может быть представлена разложением в ряд Тейлора:

Потенциальная энергия молекулы целиком имеет электромагнитную природу и обычно задается в виде суммы отдельных составляющих:

которые соответствуют следующим типам взаимодействий: и ь - потенциальная энергия валентных связей; С/„ - валентных углов; С/ ф - торсионных углов; и г - плоских групп; и пу - ван-дер-ваальсовых сил; V е1 - электростатических сил; и нь - водородных связей. Эти составляющие имеют различный функциональный вид.

Валентные связи поддерживаются за счет потенциала

где г - номер связи в молекуле; ЛГ 6 - полное число валентных связей; К Ь1 - эффективная жесткость валентной связи; г 1 - длина связи; г 0 - равновесная длина связи.

Сравнение параболического (і) и реального (2) потенциалов для валентной связи

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

Валентные углы задаются потенциалом

где I - номер валентного угла; N^ - полное число валентных углов; К,. I - эффективная упругость валентного угла; а, - значение валентного угла; а 0 , - его равновесное значение.

Энергия торсионных взаимодействий и потенциалов, соответствующих плоским группам, записывают в одинаковом виде:

где ср - номер торсионного угла; I -номер гармоники; К 0 ; - константа; # ф г - вклад гармоники в потенциал торсионного угла; л ф, - кратность гармоники. Потенциалы?/ ф и и[ отличаются константами.

Ван-дер-ваальсовы взаимодействия атомов, разделенных тремя и более валентными связями, описывают потенциалами Леннард-Джонса:

Параметры потенциала АиВ зависят от типов атомов г и у, участвующих во взаимодействии; г, = |г,- - г у |, где г, и г у - координаты взаимодействующих атомов.

Электростатические взаимодействия задают кулоновским потенциалом

где, q j - парциальные заряды на атомах; р - диэлектрическая проницаемость среды.

Водородные связи возникают и исчезают в процессе движения атомов между теми из них, которые имеют электростатические взаимодействия. Функциональный вид потенциала водородной связи похож на потенциал ван- дер-ваальсовых взаимодействий (рис. 5.5), но с более короткодействующими силами притяжения:

Водородная связь относится к специальному типу связи и обусловлена тем, что радиус иона Н + на порядок меньше, чем у других ионов. В формулах (5.39) и (5.41) имеется различие во вкладах, описывающих притяжение. Зависимость

Сравнение потенциалов для водородной связи и для ван-дер- ваальсового взаимодействия

В/г*; в (5.39) соответствует дисперсионному диполь-ди- польному взаимодействию, а величина В’/ г (Ф в (5.41) вводится, исходя из экспериментальных данных.

Отметим, что система потенциалов (5.36)-(5.41) является весьма приближенным способом задания потенциальной энергии. Его недостатки состоят в том, что энергия взаимодействия представляется в виде суммы парных сферически симметричных взаимодействий. И то и другое, вообще говоря, неверно, но это приходится пока допускать из-за отсутствия других зависимостей.

Методы поиска равновесной структуры молекулы, для которой одновременно выполняется условие равенства нулю первых производных потенциальной энергии и по всем координатам х, (ди/дх 1 = 0), могут быть разделены на две группы: минимизации с первыми производными (линейные методы) и со вторыми производными (квадратичные методы). Простые методы поиска учитывают только наклон потенциальной поверхности (т. е. первую производную), рассчитываемый численным способом или же аналитически, а более сложные методы минимизации используют как наклон, так и кривизну потенциальной поверхности (т. е. первые и вторые производные).

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

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

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

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

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

Суммирование охватывает валентные и торсионные (двугранные) углы, а также все взаимодействия между валентно-несвязанными атомами, кроме взаимодействия атомов, связанных с одним общим атомом.

Более сложные силовые поля могут включать кроме связей, указанных в уравнении (5.42), также перекрестные члены, учитывать электростатические взаимодействия и т. д. Для каждого из этих случаев имеется первое приближение, а во многих случаях разработаны и более высокие приближения. Сумма всех указанных членов называется стерической энергией молекулы.

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

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

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

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

Другой простой вариант силового поля, которое наилучшим образом соответствует принятым представлениям о природе действующих в молекуле сил, называют валентным силовым полем. Оно задается внутренними координатами, которыми обычно являются все длины связей и набор независимых валентных и торсионных (двугранных или азимутальных) углов. Это означает, что возвращающие силы действуют вдоль и поперек ковалентных связей, стремясь восстановить равновесные длины связей г, валентные углы а и торсионные углы ср. При таком выборе координат простейшее приближение, хотя и очень грубое, состоит в пренебрежении всеми недиагональными силовыми постоянными. В результате получается силовое поле, состоящее из гармонических потенциалов согласно закону Гука:

где / - функции, описывающие компоненты силового поля.

Согласно основной модели, применяемой в молекулярной механике, атомы в молекуле как бы связаны вместе отдельными независимыми пружинками, сохраняющими «естественные» величины длин связей и валентных углов. Тогда, как в случае диагонального силового поля, можно воспользоваться (по закону Гука) гармоническими потенциальными функциями, представленными уравнением (5.44) для растяжения связей и (5.45) для деформаций валентных углов:

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

имеет приемлемые свойства в определенном интервале и хорошо применима к аномально длинным связям.

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

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

• ? расстояние, при котором энергия минимальна, что связано с ван-дер-ваальсовым радиусом;
• ? глубина потенциальной ямы, что обусловлено поляризуемостью;
• ? крутизна левой части кривой, соответствующей отталкиванию, что связано с жесткостью потенциала.

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

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

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

Кулоновское взаимодействие, описываемое уравнением

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

Электростатическое диполь- дипольное взаимодействие

где?) - эффективная диэлектрическая проницаемость;

X - угол между двумя диполями р, и р у -; а, и ау - углы, образованные диполями с вектором, который их соединяет, как показано на рис. 5.6.

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

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

Разумеется, уравнение (5.50) является только первым приближением для описания торсионной энергии, выступающей в роли корректирующего члена в силовом поле. В этой косинусоидальной функции минимумы соответствуют шахматным конформациям, а максимумы - затененным. Считается, что торсионная энергия возникает вследствие отталкивания между связями, не учтенного ван-дер-ваальсовыми взаимодействиями. Альтернативным, но не полностью удовлетворительным объяснением является представление о торсионной энергии как поправке на анизотропию ван-дер-ваальсовых сил отталкивания, которое сильнее проявляется при малых, чем при больших углах к химической связи. В результате многочисленных квантовохимических работ были предложены и другие объяснения.

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

С помощью статистической термодинамики можно вычислить термически усредненную энтальпию одного кон- формера определенной молекулы, отсчитав ее от дна потенциальной ямы. Если имеется несколько конформеров, то энтальпия смеси (/Т см) находится из мольных долей (./V,) и теплот (Н:) существующих конформеров по уравнению

Мольные доли получают из распределения Больцмана

где g i - статистический вес конформера г, т. е. число идентичных конформаций; А С* - свободная энергия Гиббса без энтропийных эффектов, обусловленных наличием идентичных конформаций.

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

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

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

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

В ряде случаев удобно ввести новую величину, так называемое собственное напряжение, которое носит предварительный характер в ходе нахождения полной энергии напряжения. Собственное напряжение равно сумме показателей связи и теплоты образования наиболее устойчивого конформера, которые рассчитаны при двух следующих предположениях: 1) соединение существует в виде единственного конформера и 2) молекула не содержит открытоцепочечных связей С-С, требующих введения специальных торсионных поправок. Таким образом, полное напряжение складывается из собственного напряжения и двух дополнительных поправок. Собственное напряжение представляет интерес для одних целей, а полное - для других. Определяя эти два понятия по отдельности, можно разделить всю имеющуюся информацию, что облегчает анализ сложных систем.

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

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

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

Оптимизация геометрии в кристалле может осуществляться тремя различными способами.

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

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

Третий и наиболее мощный способ состоит в одновременной оптимизации внутри- и межмолекулярных координат.

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

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

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

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

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

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

Метод молекулярной механики успешно позволяет минимизировать энергию для больших молекулярных систем и кластеров, при разумных вычислительных затратах. Результаты расчетов методом молекулярной механики зависят от параметризации силового поля. Численное значение параметров выбирается так, чтобы получить согласие рассчитанных и экспериментальных характеристик молекулы. Параметрами являются равновесные межъядерные расстояния (длины связей) и валентные углы, а также силовые постоянные, то есть коэффициенты жесткости упругих сил, связывающих пары атомов . Простейшие модели молекулярной механики учитывают растяжения связей (Uраст.), деформацию валентных (Uдеф.) и двугранных (торсионных) углов (Uторс.), взаимодействие валентно несвязанных атомов, называемое также ван-дер-ваальсовым взаимодействием (Uвдв.), электростатические вклады (Uэл-стат.) и т.д. :

U = Uраст + Uдеф + Uторс + Uвдв + Uэл-стат (1)

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

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

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

Метод молекулярной механики неприменим: 1) для моделирования систем, свойства которых определяются электронными эффектами типа орбитальных взаимодействий и 2) в случае разрыва химических связей .

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

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

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

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

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

Другой простой вариант силового поля называется валентным силовым полем . Он наилучшим образом соответствует принятым представлениям о природе действующих в молекуле сил. Валентное силовое поле задается так называемыми внутренними координатами, а именно:

Все длины связей;

Все независимые валентные (двугранные) углы:

Все независимые торсионные (азимутальные) углы.

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

14. Достоинства и недостатки методов молекулярной механики

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

«-«Однако методы молекулярной механики могут успешно применяться лишь для сравнительно узкого класса молекулярных структур в конфигурациях, близких к равновесному состоянию.

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

2014 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Сер. 10 Вып. 1

ИНФОРМАТИКА

А. В. Райк, В. А. Клемешев

МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ МОЛЕКУЛЫ ВОДЫ С ПОВЕРХНОСТЯМИ КРИСТАЛЛОВ

Санкт-Петербургский государственный университет, 199034, Санкт-Петербург, Российская Федерация

Рассмотрена модель взаимодействия молекулы воды с кристаллической поверхностью оксида магния. Потенциальная энергия системы рассчитывалась методом молекулярной механики, учитывающим атом-атомные взаимодействия с помощью модельных полуэмпирических потенциалов. Фрагмент поверхности кристалла был представлен кластерной моделью, состоящей из конечного числа атомов, принадлежащих самой поверхности и ближайшим к ней атомным плоскостям. Были описаны различные кластерные модели, содержащие от 9 до 24 атомов оксида магния. Для них была определена равновесная геометрия системы с адсорбированной молекулой воды. Показано, что в точке равновесия молекула воды расположена на расстоянии 3 Ä от поверхности, а взаимодействие молекул воды с поверхностью носит электростатический характер. Вычислены частоты валентных колебаний адсорбированной на поверхности молекулы воды. Полученные результаты сравниваются с результатами квантово-механических расчетов. Библиогр. 15 назв. Ил. 3. Табл. 1.

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

Raik A. V., Klemeshev V. A. Modeling of interaction between a water molecule and crystal surfaces // Vestnik of St. Petersburg University. Ser. 10. Applied mathematics, computer science, control processes. 2014. Issue 1. P. 138-146.

An interaction model between a water molecule and magnesium oxide crystal surface is considered. The potential energy of the system is calculated with help the molecular mechanics method, that takes into account the atom-atom interaction using model semiempirical potentials. Fragment of the crystal surface is represented by the cluster model, that consists of a finite number of atoms belonging to the surface and the nearest atomic planes. Different cluster models containing from 9 to 24 atoms of magnesium oxide were considered. It is shown that the distance between the water molecule and the surface in equilibrium point is 3 A. It is found that the interaction is electrostatic in its nature. Stretching frequencies of the water molecule adsorbed at the surface are calculated. The results obtained are compared with the quantum-mechanical calculations. Bibliogr. 15. Il. 3. Table 1.

Keywords: molecular mechanics, intermolecular interaction potentials, adsorption, metal oxides.

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

© А. В. Райк, В. А. Клемешев, 2014

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

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

Запишем потенциальную энергию системы

Е =АЕп + ДЕ^ + ДЕШ + ДЕ7 + АЕе. (1)

Здесь ДЕд, ДЕ^, ДЕШ, ДЕ7 и ДЕе - функции, характеризующие соответственно деформацию длин химических связей, деформацию валентных углов, деформацию торсионных (двугранных) углов, невалентные взаимодействия (между химически несвязанными атомами) и электростатическое взаимодействие.

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

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

ЛЕ 1 г!2Е о

АЕЯ = Е(0) + ^ (ДАВ - ДАВ) + (ДАВ - ДАВ)2. (2)

В точке Д = До геометрия соответствует минимуму энергии, относительно которого рассчитываются все остальные энергии. Следовательно, это значение можно принять равным нулю, т. е. первый член разложения энергии в ряд Тейлора Е(0) = 0. Второй член разложения также равен нулю, так как первая производная функции в точке ее экстремума обращается в нуль. Таким образом, получаем, что потенциальная энергия зависит от третьего и высших членов разложения функции в ряд. В результате выражение (2) для энергии деформации длин связей примет вид

ДЕд = кАВ (ДДАВ)2 , (3)

где kAB - силовая постоянная, описывающая жесткость связей. Расчет энергии деформации по формуле (3) не требует больших затрат машинного времени. В случаях, когда молекулы имеют длинные связи и выходят за пределы применимости формулы (3), можно ввести дополнительный член ряда, пропорциональный (ДДАВ) .

Для описания второго слагаемого в выражении (1) - энергии деформации валентных углов - можно также использовать разложение в ряд Тейлора. Ограничимся квадратичным членом, тогда выражение для ДЕv

ДЕ^ = kABC (Д^АВ°)2 ,

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

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

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

АЕш = "у I1 + cos("-^abcd)]

При n = 1 период вращения равен 360°, при n =2 - 180° и т. д. Величины Un определяют высоту барьера вращения относительно связи C-D. В определенных молекулах при некоторых значениях n величины Un могут быть равны нулю .

Невалентные взаимодействия рассчитываются как сумма энергии отталкивания и энергии притяжения атомов. Отталкивание проявляется на малых расстояниях и может быть выражено членом в степенной зависимости от расстояния (r-12), либо в экспоненциальной форме (e-ar). Использование экспоненциальной функции приводит к тому, что полный парный потенциал взаимодействия Бакингема-Хилла стремится к -ж при r ^ 0. Кроме того, расчеты энергии молекул, включающих экспоненциальные функции, более требовательны к вычислительным ресурсам. Притяжение нейтральных атомов между собой на больших расстояниях обусловлено дисперсионными силами, главная составляющая энергии которых пропорциональна (r-6). Поэтому ДЕ7 можно записать следующим образом:

где г - расстояние между валентно-несвязанными атомами; к^ - эмпирический параметр, характеризующий интенсивность взаимодействия (глубину потенциальной ямы); коэффициент 2 - эмпирически подобранная константа.

Рис. 1. Расположение векторов при расчете электростатического потенциала, создаваемого двухатомной молекулой ЛБ в точке Р

Распределение заряда в молекуле, состоящей из электронов и ядер, неоднородно. Проанализируем вначале постоянное электрическое поле, создаваемое молекулой, поместив центр масс двухатомной молекулы АВ с зарядами ял и яв, сосредоточенными на атомах в точках -гА и +zв, в начало координат (рис. 1). Электростатический потенциал, создаваемый молекулой в точке Р, равен (ео - электрическая постоянная)

4пеЛ тл тп

_ (т2 + г"Л + 2тzA

(т2 + г"В - 2тzв

4пеот (1 + ^л/т)2 + 2(л/т)ссвв)1/2 (1 + (гв /т)2 - 2(гв /т)ссвв)1/2/

Разложив выражение (4) в степенной ряд Тейлора по обратным расстояниям при условии, что -гл ^ т и +zв ^ т, получим

(3 сся в - 1) + ...

Таким образом, в точке Р электростатический потенциал - это сумма вкладов, связанных с электростатическими моментами молекулы: я = ял + яв (полный заряд молекулы), ^ = ялzi л + Явzi в (дипольный момент) и т. д. В случае мультипольного приближения энергия электростатического взаимодействия может быть представлена в виде бесконечного ряда, включающего парные взаимодействия: заряд-заряд, заряд-диполь, диполь-диполь и др., причем все члены ряда зависят от различных степеней обратного расстояния между молекулами . Соответственно основной вклад в энергию

электростатического взаимодействия дает монополь-монопольный член (первый член разложения). Ограничившись учетом только этого взаимодействия, можно представить энергию ДЕе законом Кулона:

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

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

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

Современные программные пакеты позволяют выбрать метод, с помощью которого будет осуществляться оптимизация геометрии. Наиболее распространенными являются метод наискорейшего спуска и метод Ньютона-Рафсона . В первом случае силы, действующие на атомы, определяются методом конечных разностей, в котором для нахождения градиента в данной точке кривой сначала рассчитывают энергию в ней (рис. 2). Затем после изменения геометрии на величину бг вычисление энергии повторяется. Разность энергии бЕ в двух точках используется для нахождения градиента В методе Ньютона-Рафсона первые и вторые производные энергии находятся аналитически на основании геометрических параметров системы. Значения вторых производных, или силовых постоянных, применяются для установления положения минимума, так как они определяют кривизну потенциальной функции.

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

Рис. 2. Нахождение сил, действующих на атомы методом конечных разностей

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

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

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

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

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

Были взяты разные кластерные модели, содержащие от 9 до 324 атомов оксида магния. Для учета влияния поверхности оксида магния на адсорбированную молекулу воды рассматривались поверхностные слои размерами , , , . Число этих слоев варьировалось от 1 до 4. В качестве начальной геометрии использовался симметричный кластер с межатомным расстоянием 2 А. Выбор модели поверхностного кластерного слоя сказывался на результатах моделирования. При изучении взаимодействия одного слоя 1, 1, !, ! с молекулой воды наблюдалось сильное искажение атомов поверхности. Уже для кластерной модели двух слоев происходила заметная стабилизация геометрии, которая при увеличении слоев меняется незначительно. Атомы кислорода поверхностного слоя выступают над поверхностью (0.24 А для двух слоев и 0.17 А для четырех слоев). Очевидно, что на краях кластера искажение кристаллической структуры более заметно, поэтому молекула воды в начале оптимизации всегда помещалась над центром кластера. Расстояние от поверхности до атома кислорода молекулы воды составляло 2 А.

Рис. 3. Молекула воды над кластером М^О

В процессе оптимизации геометрии молекула воды свободно двигалась в потенциальном поле, сначала отталкиваясь от поверхности, а затем фиксируя положение равновесия. Энергетически выгодная равновесная конфигурация достигалась в точке, когда атом кислорода воды располагался над атомом магния поверхности на расстоянии 3.0-3.2 А. При увеличении модели кластера это расстояние сходится к значению 3.01 А (рис. 3). Оно хорошо совпадает с расстоянием, соответствующим точке экстремума потенциальной поверхности взаимодействия молекулы воды с кристаллом оксида магния, вычисленной методом функционала электронной плотности . Полученная в результате оптимизации геометрия расположения адсорбированной молекулы воды

над разными кластерными моделями поверхности, а именно атом кислорода воды над электроположительным катионом магния, свидетельствует о электростатическом механизме взаимодействия полярной молекулы воды с поверхностью, что хорошо передается в рамках данного силового поля. Используя такой подход, невозможно рассмотреть специфический механизм взаимодействия воды с поверхностью с возможностью образования водородных связей, как было сделано в работе в рамках метода функционала электронной плотности. Метод силового поля существенно менее требователен к вычислительным ресурсам, потому расчет подобных систем можно производить в специализированных программных пакетах даже на среднестатистическом персональном компьютере. Так, для достижения экстремума системы 4 потребовалось около 50 000 итераций и 20 мин машинного времени, для системы 2 - 12 мин. Для сравнения аналогичный квантово-механический расчет системы 2 занимает 933 ч, более 38 дней в параллельном режиме вычислений на 4 ядрах. Из данных таблицы можно видеть, каким образом происходит искажение нежесткой молекулы воды над поверхностью. Оптимизированная геометрическая структура молекулы практически не меняется. Различные кластерные модели также практически не сказываются на частоте симметричного валентного колебания молекулы воды.

Равновесная геометрия системы (А) и частота симметричного колебания Н20 (см-1)

Модель г-(Мё-Он2о) г(О-Нх) г(0-Н2) ¿(Щ-О-Нз) Ь"в

Начальная геометрия

2.00000 1.01980 1.01980 107.896

Кластер

3.49301 0.97987 0.97987 104.49 3717.45

[з х 3] 2 3.29203 0.97996 0.97994 104.52 3716.72

[ЗхЗ]4 3.22047 0.97999 0.97997 104.51 3716.44

Кластер

1 3.66293 0.97981 0.97986 104.52 3717.82

2 3.20592 0.97979 0.97985 104.49 3717.92

4 3.07547 0.97990 0.97994 104.51 3716.93

Кластер

1 3.16073 0.97996 0.97995 104.50 3716.67

2 3.07384 0.97996 0.97997 104.51 3716.54

4 3.04852 0.97997 0.97998 104.51 3716.49

Кластер

1 3.14394 0.98003 0.97993 104.49 3716.37

2 3.01869 0.97993 0.97994 104.49 3716.43

4 3.01557 0.98001 0.98002 104.52 3716.11

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

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

Литература

1. Ермаков А. И. Квантовая механика и квантовая химия: учеб. пособие. М.: Юрайт, 2010. 555 с.

2. Кларк Т. Компьютерная химия: Практ. руководство по расчетам структуры и энергии молекулы / пер. с англ. А. А. Коркина; под ред. В. С. Майстрюкова, Ю. Н. Панченко. М.: Мир, 1990. 383 с. (Clark T. A Handbook of computational chemistry.)

3. Цирельсон В. Г. Квантовая химия. Молекулы, молекулярные системы и твердые тела. М.: БИНОМ, 2010. 496 с.

4. Каплан И. Г. Межмолекулярные взаимодействия. Физическая интерпретация, компьютерные расчеты и модельные потенциалы. М.: БИНОМ, 2012. 394 с.

5. Каплан И. Г. Введение в теорию межмолекулярных взаимодействий. М.: Наука, Гл. ред. физ.-мат. лит., 1982. 312 с.

6. Rogers D. Computational chemistry using the PC. Ed. 3. Hoboken, New Jersey: John Wiley & Sons, 2003. 363 p.

7. Бедрина М. Е., Егоров Н. В., Клемешев В. А. Моделирование наноструктур на высокопроизводительном вычислительном комплексе // Вестн. С.-Петерб. ун-та. Сер. 4: Физика, химия. 2010. Вып. 4. C. 136-140.

8. Stewart J. J. P. Optimization of parameters for semiempirical methods. I. Method // Journal of Computational Chemistry. 1989. Vol. 10, N 2. P. 209-220.

9. Bandura A. V., Sykes D. G. Adsorption of Water on the TiO2 (Rutile) (110) Surface: A Comparison of Periodic and Embedded Cluster Calculations // J. Phys. Chem. B 2004. Vol. 108, N 23. P. 7844-7853.

10. Fox H., Gillan M. J., Horsfield A. P. Methods for calculating the desorption rate of an isolated molecule from a surface: Water on Mg0(001) // Surface Science. 2007. Vol. 601. P. 5016-5025.

11. Bellert D, Burns K. L., Wampler R. e. a. An accurate determination of the ionization energy of the MgO molecule // Chem. Phys. Lett. 2000. Vol. 322. P. 41-44.

12. Абаренков И. В., Третьяк В. М., Тулуб А. В. О механизме диссоциации молекул воды на (100) поверхности кристалла MgO в модели парных неэмпирических потенциалов // Хим. физика. 1985. Т. 4, № 4. С. 974-980.

13. Егоров Н. В., Денисов В. П., Ермошина М. С. Моделирование взаимодействия заряженной частицы с поверхностью электрода сложной конфигурации // Поверхность. 2005. № 7. C. 60-63.

14. Райк А. В., Егоров Н. В., Бедрина М. Е. Моделирование потенциалов межмолекулярного взаимодействия // Вестн. С.-Петерб. ун-та. Cер. 10: Прикладная математика, информатика, процессы управления. 2012. Вып. 3. С. 79-87.

15. Райк А. В., Бедрина М. Е. Моделирование процесса адсорбции воды на поверхности кристаллов // Вестн. С.-Петерб. ун-та. Cер. 10: Прикладная математика, информатика, процессы управления. 2011. Вып. 2. С. 67-76.

Контактная информация

Райк Алексей Владимирович - аспирант; e-mail: [email protected]

Клемешев Владимир Алексеевич - кандидат физико-математических наук, доцент; e-mail: v. klemeshev@spbu. ru

Raik Aleksey Vladimirovich - post-graduent student, St. Petersburg State University, 199034, St. Petersburg, Russia Federation; e-mail: [email protected]

Klemeshev Vladimir Alekseevich - candidate of physical and mathematical sciences, associate professor, St. Petersburg State University, 199034, St. Petersburg, Russia Federation; e-mail: [email protected]