Математическое моделирование эффектов сольватации. Разработка модели и комплекса программ тема автореферата и диссертации по химии, 02.00.04 ВАК РФ

НАЦИОНАЛЬНАЯ АКАДЕМИЯ НАУК РЕСПУБЛИКИ КАЗАХСТАН

ИНСТИТУТ ХИМИЧЕСКИХ НАУК им. А.Б.БЕКТУРОВА

На правах рукописи УДК 519.688:(541.8+539.196)

РУБАНЮК Надейгда Николаевна

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЭФФЕКТОВ

СОЛЬВАТАЦИИ. РАЗРАБОТКА МОДЕЛИ И КОМПЛЕКСА ПРОГРАММ.

02.00.04 - физическая химия

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата химических наук

АЛ МАТЫ,1995

Работа выполнена в Институте химических наук им.А.Б.Бектурова Национальной академии наук Республики Казахстан.

Научный руководитель - кандидат химических наук, с.н.с.

Габдракипов В.З.

Официальные оппоненты - доктор химических наук,

профессор Оспанов Х.К.

кандидат химических наук, доцент Дильмухамбетов Е.Е.

Ведущая организация - Институт органического катализа и электрохимии им.Д.В.Сокольского HAH PK, г.Алматы.

Защита состоится " 72." января 1996 г. в 10 часов на заседании Специализированного совета К 53. ¡8.02 при Институте химических наук 1Ш.А.Б.Бектурова HAH PK по адресу: 480100, г.Алматы, ул.Вапиханова, 106.

С диссертацией можно ознакомиться в библиотеке ИХН им.А.Б.Бектурова HAH PK.

Автореферат разослан " " декабря 1995 г.

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

Б.С.Балапанова

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ.

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

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

Связь с планом основных научных работ. Диссертационная работа выполнена в соответствии с темой научно-исследовательской работы Института химических наук им. А.Б.Бектурова "Создать научные основы новых направлений переработки фосфатного сырья на продукты химизации сельского хозяйства и неорганические материалы" (К госрегистрации 01890052171).

Научная новизна работы. В настоящей работе предлагается простой, но достаточно строгий и не требующий огромных затрат ма-

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

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

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

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

Автор выносит на защиту:

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

- комплекс программ, реализующий модель учета эффектов сольватации;

- потенциал водородной связи;

- результаты моделирования сольватных оболочек однозарядных ионов.

Публикация и аппробация работы. По материалам работы опубликовано 6 работ. Результаты работы докладывались на VII Всесоюзной конференции "Использование вычислительных машин в химических исследованиях и спектроскопии молекул" (Рига ,1986 г.), II Республиканской конференции "Проблемы вычислительной математики и автоматизации научных исследований" (Алма-Ата, 1988 г.), VI Всесоюзной конференции молодых ученых по физической химии "Физхимия -90" (Москва, 1990), Всесоюзном координационном совещании по квантовой химии (Новосибирск, 1990).

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

рисунков. Список литературы включает 180 источников.

ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ.

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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И АЛГОРИТМ.

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

В диэлектрическом континиуме с макроскопической диэлектрической постоянной £о вырезается сферическая полость с е0 =1. В центр полости помещается молекула (ион) растворенного вещества, вокруг которой в пределах полости размещается конечное число молекул растворителя (до 100). Радиус полости определяется как расстояние от

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

Межмолекулярные взаимодействия вычисляются в аддитивном (парном) приближении и описываются методом атом-атомных потенциалов, включая кулоновские:

= (1)

I }

= (2)

и

где IIу - атом-атомный потенциал Леннард-Джонса, Букингема или Смитта-Теккера, (7ЭЛ - кулоновская часть потенциала (2), г у - расстояние между атомами г и заряды на атомах 1 и Кулоновские взаимодействия на малых расстояниях, когда волновые функции молекул начинают сколь-нибудь существенно перекрываться, корректируются функцией переключателя Б (г ¡¡), записанного в форме кубической параболы (переключатель Стиллинжера), квартичной параболы или синусоидального переключателя.

Взаимодействие атома любой молекулы внутри полости со средой (диэлектрическим конпшиумом) учитывается в приближении теории реактивного поля по формуле Кирквуда (3), что эквивалентно методу зарядов изображения, либо по формуле Онзагера-Бетчера (4) (метод наведенного диполя)

1 А Я

<7/ +

г-г, | а ^ е(/ + 1) + /

s-1 (M-tiJ

2e + 1 o3

где энергия взаимодействия молекулы А с континиумом, г, и г, -радиус-вектора атомов / и _/, д - радиус полости, // - дипольный момент молекулы А, М - векторная сумма дипольных моментов всех молекул, находящихся в полости, Рг- полиномы Лежандра степени I.

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

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

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

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

r\\P~S(r)rw (5)

где rw/1 - принятый радиус атома водорода, S(r)~ функция, ана-

(5)

Etot - UdB + U,

рп

(6)

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

Описанный алгоритм реализован в виде комцлекса программ SOLUT на IBM PC/AT.

ОПИСАНИЕ ПРОГРАММНОГО КОМПЛЕКСА SOLUT.

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

ственные ограничения: максимальное число молекул в системе - 100, максимальное число атомов в одной молекуле - 100, максимальное число типов молекул в системе - 10. Все программы пакета БОЬиТ написании на языке Фортран 77 с последующим компилированием этих программ компилятором фирмы МкгоБоЛ версий 5.0 и выше. Размер загрузочного модуля комплекса БОШТ составляет 310 кБ.

ПАРАМЕТРИЗАЦИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ.

Для проверки корректности используемых приближений и оценки работоспособности предложенной дискретно-континуальной модели сольватной оболочки было проведено моделирование следующих систем: (Н20)п, М+ (НгО)п, А- (Н20)п, где п=1,2...,8,35; М+ - однозарядные ионы щелочных металлов (1л+, К+, Сб+); А- - однозарядные ионы галоидов (Р-, С1 , Вг, 1-).Всего 90 систем.

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

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

(3) и через наведенный диполь (4) будет одинаковой либо при условии, что все гибридные диполи моментов равны нулю, либо когда заряды на атомах переопределены так,чтобы включать в себя гибридные диполи.

Заряды в молекуле воды с целью включения в зарядовый диполь гибридных диполей переопределяли путем оптимизации функционала дипольного момента методом Давидона-Флетчера-Пауэла (H20-3q).

При расчете электростатических взаимодействий между молекулами необходимо учитывать энергию взаимодействия между остаточными зарядами атомов и гибридными диполями, а также между самими гибридными диполями. Для этого каждый гибридный диполь можно заменить соответствующими зарядами, помещая их на рассматриваемом атоме, или на соседних с ним атомах, или на псевдоатомах. Физический смысл этого ассимметричного заряда - моделирование направленности неподеленной электронной пары (НЭП). Мы выбирали НЭП таким образом, чтобы воспроизводить полный дипольный момент молекулы, поэтому длину n-пары и ее заряд мы калибровали по дипольному моменту. Для нахождения НЭП использовали либо локализованные молекулярные орбитали (J1MO), либо валентные атомные орбитали (BAO). Преобразуя базис АО к базису валентно-активных орбиталей (BAO) и далее к базису оптимальных гибридных орбиталей (ОГО), ориентированных по неподеленным парам электронов и связям в молекуле, определяли угол между двумя ОГО, относящимися к НЭП. В случае Л МО канонические молекулярные орбитали преобразовали в локализованные МО (ЛМО), затем в ЛМО находили те, которые описывают НЭП и коэффициенты ЛМО при гетероатоме использовали для определения направления НЭП. Длину вектора НЭП определяли через среднее квантовомеханическое <г2>,/2 , заряды на НЭП определяли из условия электронейтральности суммарного заряда в молекуле fH20-5q).Геометрия и зарядовое распределение для обеих моделей приведено в табл.1.

Таблица 1.

Геометрия (Я нм, а, ц град.) и заряды ((¡>, ед.электрона) в трехточечной (Н20-Зд) и пятиточечной (Н20-5д) молекулах воды

Яо-н Яо-„ ССн-о-н <5о <3н

НаО-^ 0.961 - 104 - -0.660 0.330 -

Н20-5Я 0.961 0.875 104 127.8 -0.030 0.191 -0.176

МОДЕЛИРОВАНИЕ ДИМЕЮВ И МАЛЫХ КЛАСТЕРОВ ВОДЫ.

Моделирование структуры и энергетики димера и малых кластеров воды в программном комплексе 80ШТ проводилось для двух моделей воды-трехточечной и пятиточечной. Полученная в результате оптимизации геометрия димера является оптимальной и соответствует одному из локальных минимумов на потенциальной поверхности, описывающей взаимодействие двух молекул воды. Различные геометрические формы димера удается получить, меняя стартовую точку оптимизации. Отличие в энергии между различными конфигурациями димера невелико и составляет ±0.5 ккал/моль. Наиболее глубокие минимумы энергии на потенциальной поверхности соответствуют трем известным формам димера - линейной Е=-5.03 (-5.59) ккал/моль, циклической Е=-4.60 (-5.16) ккал/моль и бифуркатной Е—3.65 (-3.77) ккал/моль (в скобках значения для пятиточечной модели молекулы воды). Для обеих моделей воды самой компактной является циклическая структура: 11о...о=0.246 (0.256) нм .наибольшее расстояние Я0...о соответст вует бифуркатному димеру-0.291 (0.293) нм,в линейном да-мере Яо..о =0.281 (0.280) нм, что несколько меньше экспериментального.

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

а более сложные системы, свойства соединений заметно изменяются. Нами проведено моделирование малых ассоциатов воды, включающих до 8 молекул. Оптимальные структуры ассоциатов были получены путем минимизации полной энергии взаимодействия всех молекул, входящих в кластер. Как и в случае димера воды, на поверхности потенциальной энергии каждого ассоциата имеется большое количество близких по глубине потенциальных ям, соответствующих точкам локальных минимумов, т.о. равновесные структуры водных кластеров не являются строго фиксированными из-за существования близких по энергиям конфигураций, отвечающих одному и тому же ассоциату. Но среди всех минимумов энергий наиболее глубокие соответствуют уплощенным кольцевым структурам, в которых каждая молекула воды выступает одновременно и донором, и акцептором в Н-связи. Образование более крупных кластеров сопровождается большим энергетическим выигрышем в расчете на одну молекулу благодаря возникновению замкнутых конфигураций с Н-связями. Полученные нами в результате моделирования структуры не обладают симметрией и являются нерегулярными. С увеличением размера кластеров наблюдается возрастание их стабильности, одновременно увеличивается и энергия, приходящаяся на одну молекулу НгО. Среднее расстояние Я0...0 с укрупнением ассоциатов (для п>2) практически не меняется и по величине ближе к межкислородному расстоянию в циклическом димере. В термодинамическом смысле наличие чуть более глубокого минимума потенциальной энергии у одной структуры не свидетельствует о ее более вероятном образовании из мономеров, и решение вопроса об устойчивости ассоциата должно проводиться с учетом энергетики и энтропийного фактора. Анализ водородных связей в малых кластерах воды показал, что, как правило, каждая молекула воды стремится ориентироваться таким образом, чтобы образовать максимальное количество К-связей с соседями, при этом одновременно может являться как донором, так и акцептором протона. Н-связи, образованные одной

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

РАСЧЕТЫ КЛАСТЕРОВ КАТИОНОВ ЩЕЛОЧНЫХ МЕТАЛЛОВ И

ГАЛОИД АНИОНОВ.

ОПРЕДЕЛЕНИЕ СТРУКТУРНЫХ И ЭНЕРГЕТИЧЕСКИХ

ХАРАКТЕРИСТИК СОЛЬВАТНЫХ ОБОЛОЧЕК ИОНОВ.

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

При моделировании гидратации малых кластеров однозарядных ионов, включающих до 8 молекул НгО, либо первоначально проводилось "намораживание" очередной молекулы ШО на уже существующий комплекс и далее осуществлялась процедура полной оптимизации энергии созданной "капли", либо все молекулы Н2О (п=1,8) размеща-

лись вокруг иона и путем полной оптимизации энергии определяли наиболее оптимальную структуру "капли". В ходе оптимизации для всех кластеров было найдено большое количество структур, соответствующих локальным минимумам энергии, которые не существенно отличаются друг от друга как энергетическими, так и структурными характеристиками. Сравнение результатов нашего расчета с экспериментальными и расчетными данными показывает хорошую корреляцию для всех ионов. Для всех однозарядных ионов экспериментальные значения ЛНп-1,п уменьшаются с ростом Пшо, а при постоянном Пню изменение энтальпии максимально в случае Тг иИ-и минимально для Сб+ и 1-, т.е. коррелирует с ионным радиусом. Расчитанные нами энергии последовательной гидратации сохраняют ту же тенденцию изменения с ростом числа молекул НгО.

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

Для определения структурных и энергетических характеристик

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

АН к Е (п+юп) - Е (п) или АН я Е (п+хоп) - Е (п-Н), (7)

где п - число молекул воды в кластере. Поиск оптимальных конфигураций моделируемых систем осуществлялся процедурой минимизации полной энергии системы. При этом, как и в случае малых кластеров, было обнаружено большое количество структур, близких по энергии. В процессе моделирования поиск водородных связей в сольватных оболочках ионов осуществлялся по геометрическим (Яо._о ^0.32 нм, Ы0...н< Ко...о, а (0-Н...0)<90°) и энергетическому критерию (Е $ -1 ккал/моль), при этом определялся тип координации молекул и проводился анализ геометрических и энергетических характеристик Н-связей. Результаты моделирования сольватных оболочек однозарядных ионов и воды приведены в табл.2 и 3. Совпадение расчитан-ных с экспериментальными данными хорошее: в среднем отклонения вычисленных нами энергий гидратации от экспериментально найденных АН составляют ±5.1 ккал/моль, предложенная нами математическая модель позволяет вычислять энергии гидратации с точностью до 94.5%. Определенные нами радиусы сольватных оболочек и координационные

Таблица 2.

Энергии (Е) и энтальпии (эксп.,ДН) гидратации ионов, вклад электростатики (Ея) и реактивного поля (Ерп) (% от Е), средние энергии взаимодействия молекул Н2О 1-ой сольватной оболочки с ионом СЕ:) радиусы кавитации (Яс), радиусы первых сольватных оболочек (Яг..о) и координационные числа (£) ионов. Энергии в ккал/моль, радиусы в им. В скобках - литературные данные.

-Е Ее Ерп -ДН 11с -Е1 Ъ

Н20 - 93.5 3. - 1.46 - - -

135.3 92.2 4.3 128.4 1.29 31.7 0.182(0.190-0.225) 5(4-6)

100.6 90.2 4.8 101. 1.28 22.0 0.223 (0.229-0.250) 6(5-7)

К+ 90.1 90.5 4.7 80.9 1.34 16.1 0.261 (0.260-0.295) 7(5-7)

11Ъ+ 75.8 89.6 5.2 75.1 1.29 14.5 0.279 (~ 0.29) 6(6-14)

Ол+ 73.6 80.2 6.0 69. 1.18 12.0 0.301 (0.295-0.322) 7 (7-8)

и- 113.1 89.5 5.1 120.1 1.17 17.5 0.261 (0.260-0.322) 7(4-8)

С1 97.0 88.6 5.4 86.3 1.19 14.2 0.301 (0.270-0.325) 7-9 (6-8)

Вт- 76.0 87.5 5.8 79. 1.20 10.6 0.325(0.312-0.345) 10 (4-9)

У 74.6 89.4 4.4 69.2 1.29 9.3 0.344 (0.300-0.376) 10 (6-9)

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

Таблица 3.

Общее число Н-связей в кластере (Кн), распределение их по типам (Ь-линейные, В-бифуркатные, С-циклические) и феднее число Н-связей на молекулу воды (<№1>), средняя длина (Но«, нм) и средняя энергия Н-связи (Ен, ккал/моль) в гидратных кластерах ионов.

ш Ь В С <№1> К.О-0 Ен

ц+ 76 56 8 12 2.24 0.273 -3.8

80 24 32 ' 24 2.35 0.268 -4.5

К+ 100 28 28 44 2.94 0.269 -4.0

яъ+ 106 28 32 46 3.12 0.264 -4.0

102 38 28 36 3.29 0.270 -3.7

Р 90 34 16 40 2.65 0.270 -3.6

С1- 100 22 24 54 3.24 0.266 -3.6

Вг- 104 24 24 56 3.06 0.269 -3.6

л- 118 18 48 52 3.47 0.265 -3.7

шо 124 36 40 48 3.54 0.270 -4.0

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

Результаты анализа кластеров на водородные связи приведены в табл.3. В процессе анализа подсчитывались частоты встречаемости молекул Н2О, образующих разное количество водородных связей и определялся тип координации. Всего обнаружено 14 различных типов координации молекул, обозначаемых Дпат, где п - число водородных связей, в которых молекула воды участвует в качестве донора, т - в качестве акцептора протонов.С увеличением ионного радиуса в рядах Хл-Сб и Р-1 увеличивается общее число водородных связей и соответственно среднее число Н-связей, приходящихся на одну молекулу воды в кластере,уменьшается относительное число линейных Н-связей, при этом средняя энергия и средняя длина Н-связи остаются практически неизменными.

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

ВЫВОДЫ.

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

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

2. Математическая модель реализована в форме программного комплекса ЗОЫЛ. Общее количество программ и подпрограмм в комплексе-176. Комплекс снабжен 1рафическим интерфейсом, отображающим на экране процесс образования и оптимизации сольватной оболочки.

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

4. Предложен простой потенциал для описания водородной связи в виде зависимости ван-дер-ваальсова радиуса атома водорода от расстояния Х...Н (Х=М, О, И, С1, Вг, I). С этим потенциалом хорошо воспроизводятся экспериментальные энергии Н-связей, расстояния О...О и геометрия комплексов с водородной связью.

5. Изучены эффекты сольватации однозарядных ионов 1д+, Ыа+, К+, ЯЬ+, Сз+, Р , С1-, Вг, I- в воде. Получено хорошее согласие с экспериментом расчитанных энтальпий гидратации, радиусов сольватных оболочек и координационных чисел. Анализ водородных связей в изученных системах хорошо согласуется. с результатами, полученными методами молекулярной динамики и Монте-Карло.

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

точно проста и может быть применима к широкому ряду задач физической, органической, неорганической и биоорганической химии.

1. Рубанкж H.H., Габдракипов В.З. Теоретическое изучение эффектов сольватации.//Изв. АН КазССР,сер.хим.-1990.-М 6.-С. 16-20.

2. Рубанюк H.H., Габдракипов В.З. Математическая модель соль-ватной оболочки.Комплекс программ.//Изв.АН КазССР,сер.хим.-1990. - N 6.- С.20-22.

3. Рубанюк H.H., Шпнльберг И.Г., Габдракипов В.З. Программный комплекс для расчета координат атомов и некоторых физических характеристик молекулы,получения аксонометрической проекции на графопостроителе с помощью мини-ЭВМ СМ-3. //Использование вычислительных машин в химических исследованиях и спектроскопии молекул: Тез.докл.VII Всесоюзн.конф.21-23 ноября 1986.-Рига,1986,-

4. Рубанюк H.H., Габдракипов В.З. Программный комплекс для расчета координат атомов и некоторых физических характеристик молекул. // Проблемы вычислительной математики и автоматизации научных исследований: Тез.докл.И Республ. конф. 2-6 ноября 1988 г. -Алма-Ата: Наука, 1988. - Т.1.- С.88.

5. Габдракипов В.З.,Рубанюк Н.Н.Пакет подрограмм по аналитической геометрии и векторной алгебре. // Проблемы вычислительной математики и автоматизации научных исследовашш: Тез.докл.И Рес-публ.конф. 2-6 ноября 1988 г.-Алма-Ата: Наука,1988. - Т.1.- С.ЗЗ.

6. Рубанюк H.H.,Габдракипов В.З. Квантово-химическое обеспечение мини-ЭВМ Tima СМ-3 и СМ-4.Контроль расчета геометрии молекул и валентно-активные орбитали.// Изв.АН КазССР, сер.хим.,1991. -N 1.-С.23-24.

ПО ТЕМЕ ДИССЕРТАЦИИ ОПУБЛИКОВАНЫ СЛЕДУЮЩИЕ РАБОТЫ:

С.241.

РЕЗЮМЕ

Рубанюк Надежда Николайкызы СОЛЬВАТ ЭФФЕКТЮ1Н МАТЕМАТИКАЛЫК; МОДЕЛЬДЕУ ПРОГРАММА ЖИЫНЫН ЖвНЕ МОДЕЛ1Н ЖАСАУ

Химия гылыиыпыц кандидат дэрежесше 1здену 02.00.04 - физикалык, химия

Бул диссертацияда сольват зффекттстнщ ертлгеи молекула-лармен иондардыц кдеиетше осср{ тексершген. К,ойылгап мадсатты шешу угаш, жай, бхрак, ете далме—дал, машиналык, уак,ытты коп к,ажет етнейтш, сольват зффектасш есепке алатын математикалык, дискрет-континуальды моделш курайгын эдас усынылды. Сольват кдбыгын куру жэне сольват жиыныныд туракты кодфигурациясын табу оптимизация одштерг аркылы жургазйда. Сутектш байланыстарды суроттеу ушш модифика-цияланган атом-атомдьщ потенциал усынылды. Математикалык модель программа жиыны туронде жузеге асьгрылды. Ол сольват г;абыгыныц кещстштеп жэне энергетикальщ курылысын тексе-руге, модельдеугс, сонымеп б!рге молекулалардыц немесе реа-генттердш, потенциялдык, бетшде, ортамен эрекеттесулертнщ асершен пайда болатьш энергияны тузетуда есептеуге, сольват жиыпындагы сутектл байланыстардыц курылысымен энергиясын талдауга мумкшдш бередь Табылган модель бойынша, б1р за-рядты иондардыц Ы+, N8+, К+, М>+, Сз+, Б-, СЬ, Вг-, Л судагы сольват эффектиюр1 тексерьпда. Есептелген гидратацияныц энталь-пиясымен координация сандары экспериментпен жак,сы уйлестт. Тексершген жуйедег! сутектш байланыстарды талдау молекула-лык, динамика жэне Монте-Карло одастер! аркылы алынган натажелерге сэйкес келда. ¥сыпылып отырган математикалык, модель/ц физикалык, органикалык,, бейорганикалык, жэне био-химиялык, мэселелердо шешу ушш крлдануга болады.

RESUME

Rubanyuk Nadezhda Nikolaevna MATHEMATIC MODELLING OF SOLVATATION EFFECTS.

DEVELOPMENT OF THE MODEL AND COMPLEX OF PROGRAMMS Dissertation on graduation of candidate of chemical sciences 02.00.04 - physical chemistry

This work is dedicated to the study of solvatation effects upon the properties of the dissolved molecules and ions. A simple but rather strict method of building a mathematical discrete-continuous model of accounting for the solvatation effects, not demanding large expenses of machine time has been proposed. The building of solvatation envelope and the finding of stable configuration of the solvatation complex are fulfilled by optimitaion methods. For the description of hydrogen bonds a modified atom-atom potential is suggested. The mathematic model is realised in the form of the programm complex and it allows to model and study the spatial and energetic strusture of solvation envelopes, to calculate adjustments to energies on the potential surface of the molecule or reagents, stipulated by the interaction with the media, to analise the structure and energetics of hydrogen bonds in the solvatation complexes. With the help of the developed model there has been conducted the study of solvatation effects of ont-charge ions Li+, Na+, K+, Rb+, Cs+, F-, CI-, Br-, J- in water. A good correspondence with the experiment of calculated enthalpies of hydratation and coordination numbers is obtained. An analysis of hydrogen bonds in the studied systems concords well with the results, obtained by method of molecular dynamics and Monte-Carlo methods.The suggested mathematical model can be applied to a wide range of tasks of physical, organic, inorganic and bioorganic chemistry.

Подписано в печать 5.12.1995 г. Формат 60 х 84 1/i6- Бум. тип №1. Усл. п.л. 1,5.

Тираж 100 экз. Заказ № 62. Издано в Компьютерном Центре

Института Теоретической и Прикладной Математики HAH РК 480021, Алматы, ул. Пушкина, 125