Содержание к диссертации
Введение
Глава 1. Обзор литературы 10
1.1. Классические силовые поля 10
1.1.1. Стандартный вид невалентных взаимодействий 10
1.1.2. Валентные взаимодействия 12
1.1.3. Явные и неявные модели воды 13
1.2. Поляризуемые силовые поля 14
1.2.1. Модели с точечными диполями 14
1.2.2. Модель флуктуирующих зарядов 15
1.2.3. Модель на основе осцилляторов Друде 16
1.2.4. Поляризуемость на основе модели Толе 16
1.2.5. Возможные улучшения поляризуемой модели 18
1.3. Методы получения параметров СП из экспериментальных данных.. 19
1.3.1. Методы моделирования различных агрегатных состояний вещества 21
1.4. Квантовые вычисления и их использование для построения СП 23
1.4.1. Уравнение Шредингера 23
1.4.2. Метод Хартри-Фока 24
1.4.3. Наборы базисных функций 26
1.4.4. Метод МР2 и его модификации 27
1.4.5. Другие квантовые методы и их применение 29
1.5. Функциональная форма СП QMPFF3 32
1.5.1. Зарядовая плотность QMPFF3 32
1.5.2. Электростатическое взаимодействие 36
1.5.3. Обменно-отталкивающее взаимодействие 37
1.5.4. Дисперсионное взаимодействие 37
1.5.5. Поляризационное взаимодействие 38
1.6. Параметризация QMPFF3 при помощи квантово-механических данных 38
1.6.1. Квантовые вычисления для построения QMPFF3 38
1.6.2. Декомпозиция МР2 энергии для выделения четырех компонент энергии 40
1.7. Типизация и процедура параметризации QMPFF3 41
Глава 2. Создание силового поля FFSol, неявно учитывающего водное окружение молекул 46
2.1 Поле FFSol и необходимые для его параметризации экспериментальные данные 46
2.1.1 Функциональная форма силового поля FFSol 46
2.1.2 Свободная энергия растворения молекул кристалла в воде 48
2.1.3 Вклад заторможенных в кристалле степеней свободы 49
2.1.4 Необходимые экспериментальные данные 50
2.1.5 Создание базы данных CRAFT 51
2.2 Оптимизация энергии кристалла 54
2.2.1 Типизация и заряды СП 54
2.2.2 Энергия кристаллической ячейки и ее минимизация 55
2.3 Получение невалентных параметров СП 59
2.3.1 Функция невязки и ее минимизация 59
2.3.2 Процедура получения параметров 61
2.4 Тестирование силового поля 63
2.4.1 Параметры вспомогательного силового поля FFSubl и тестирование их качества 63
2.4.2 Параметры «водного» силового поля FFSol и тестирование их качества 65
Глава 3. Уточнение параметров взаимодействий ароматического углерода в СПQMPFF3 70
3.1 Коррекция параметров ароматического углерода в СП QMPFF3 72
3.1.1 Предварительная параметризация на основе МР2 метода 72
3.1.2 Коррекция дисперсионного параметра с помощью данных CCSD(T) 78
3.2 Подтверждение качества коррекции 80
3.2.1 Газовая фаза 80
3.2.2 Жидкая фаза 82
3.2.3 Кристаллическая фаза 85
3.3 Качество силового поля определяется качеством лежащих в его основе квантово-механических данных 88
Глава 4. Тестирование силового поля QMPFF3 90
4.1 Предсказание QMPFF3 свойств молекул 90
4.1.1 Поляризуемость ' 92
4.1.2 Дипольный момент 93
4.2 Энергия Гиббса димеризации в газах 94
4.3 Моделирование однородных жидкостей и сравнение с другими СП 96
4.4 Оптимизация кристаллов и сравнение с другими СП 99
4.5 Вакуумная минимизация белков как тест силового поля 102
Заключение 108
Выводы 111
Благодарности 112
Приложение
- Поляризуемые силовые поля
- Функциональная форма СП QMPFF3
- Оптимизация энергии кристалла
- Параметры вспомогательного силового поля FFSubl и тестирование их качества
Введение к работе
В настоящее время сфера биотехнологий и создания новых материалов приближается к новому технологическому прорыву. Для предсказания новых свойств материалов, для моделирования молекулярных, и, в частности, биомак-ромолекулярных систем все чаще обращаются к компьютерному эксперименту, использующему силовые поля. Силовым полем (СП) называют набор потенциальных функций и параметров взаимодействий, описывающих валентное и невалентное взаимодействие атомов в молекулах и между молекулами [1]. С помощью молекулярной динамики (МД) можно моделировать как простые жидкости [2-7], так и сложные системы, включающие несколько макромолекул и растворитель [8]. В последнем случае, реалистичность описания требует учета очень большого числа молекул растворителя, что - если молекулы рассматриваются в явном виде - резко, на порядки увеличивает время расчетов.
Таким образом, разработка силовых полей наталкивается на две принципиальные трудности. С одной стороны, СП должно работать достаточно быстро. С другой стороны, оно должно достаточно точно описывать рассматриваемые взаимодействия. Точность описания разнообразных систем с помощью МД зависит от качества СП. Высокая точность описания систем необходима, в особенности, в том случае, когда нужно выделить одну лучшую структуру из большого множества: сотню лучших лигандов из миллиона, потенциально способных связаться с целевым белком [9-12]; стабильную нативную структуру белка среди множества других [13]; наиболее подходящие, с точки зрения протонной проводимости, сульфат-ионные (-R-SO3) боковые группы в перфторполимерной мембране [14, 15] и т.п.
Точность современных СП, даже модифицированных специально под белки, до сих пор не позволяет описать моделируемые белки в воде при помощи МД с расхождением (по Са атомам) от первоначальной, экспериментально-определенной структуры менее, чем на 1.5 А [16]. Последний раунд соревнования CASP по предсказанию белковых структур по их аминокислотным последовательностям не показал существенного прогресса за последние годы, следова- тельно, такие предсказания все еще нуждаются в более точных, чем существующие, силовых полях [17]. Последнее десятилетие было посвящено разработке поляризуемых СП, в которых учитывается влияние электрического поля на поляризуемые атомные облака, так как только в таких системах может возникнуть корректная диэлектрическая проницаемость [18], а также другие эффекты, не являющиеся попарно-аддитивными [19-21], что является необходимым условием для моделирования таких сложных гетерогенных объектов, как белки.
Недавно было создано поляризуемое СП QMPFF3 общего назначения (т.е. для широкого круга задач), которое базируется исключительно на квантовых данных, что позволяет ему не зависеть от все еще недостаточного количества экспериментальных данных. QMPFF3 имеет сложную, физически-обоснованную функциональную форму. Чем сложнее и физически-обоснованнее функциональная форма силового поля, тем, потенциально, точнее может оно описать энергию взаимодействия в различных системах, и тем больше требований предъявляется к точности подбора численных параметров такого СП.
В частности, силовое поле QMPFF3 не позволяет достоверно описывать ароматические структуры, а также полиароматические углеводороды (ПАУ) [22]. Правильное описание взаимодействия ПАУ является существенной частью разработки СП общего назначения, так как ПАУ имеют большое значение для решения широкого спектра научных и технологических задач (например, разработка лекарств [12], самообразование графитных нанопроволок [23] и т.п.) и являются частями белков, ДНК, РНК, многих лекарств, наноматериалов и т.п. Большую важность ПАУ представляют и для медицинских и экологических наук из-за своей существенной канцерогенности и токсичности [24].
При исправлении недостатка, связанного с плохим описанием ароматических структур, и после качественного тестирования, которое сможет подтвердить хорошую предсказательную силу СП QMPFF3, в том числе и для белков, его можно было бы использовать для задач моделирования макромолекул, а также для предсказания свободных энергий связывания комплексов белок-лиганд [12].
Детальные и сложные СП обычно требуют существенно больше времени для вычислений, чем простые. Поэтому для МД моделирования сворачивания белков (притом, что, характерное время их сворачивания - секунды-миллисекунды, а моделирование таких процессов занимает многие тысячи часов) необходимы «быстрые» модели, желательно учитывающие воду неявно, и при этом правильно описывающие все основные компоненты взаимодействий. Несколько моделей взаимодействий с неявным учетом воды уже было представлено до сегодняшнего дня [25-29], однако до сих пор они не столь точны, чтобы использовать их для выделения нативной структуры белка. Поэтому очень интересной и важной задачей представляется создание такого силового поля для невалентных взаимодействий, которое систематически учитывало бы взаимодействие с водой в неявном виде (так, как учитывается диэлектрическая проницаемость в макроскопических моделях веществ), т.е. было бы основано на принципиально новом для атомных силовых полей подходе.
Целью данной работы является построение силовых полей для оценки взаимодействия молекул в вакууме и в неявно учитываемой воде для моделирования белков и их взаимодействий с другими молекулами.
Соответственно, были поставлены следующие задачи: 1) создать базу кристаллографических и термодинамических данных для определения параметров взаимодействия молекул в воде; 2) построить, на основании сведений из данной базы, силовое поле для оценки взаимодействий в неявно заданном водном окружении (назовем его FFSol); 3) улучшить параметризацию ароматических молекул в существующем поляризуемом силовом поле QMPFF3, которое базируется на квантово-механических данных, и провести тестирование результирующего поля.
Поляризуемые силовые поля
Модели с точечными вращающимися диполями (предполагается что своим вращением они реагируют на внешнее электрическое поле) на атомных центрах с поляризуемостью а, которые были предложены ранее всего, позволяют рассчитывать поляризуемые системы достаточно просто, но приводят к так называемой поляризационной катастрофе, когда, при сближении двух диполей с поляризуе / \і/б мостями а, и а, на расстояние rcntlcal =\4аха}) , эти два диполя-падают друг на друга [1, 81], что приводит к нестабильностям при расчете МД. Если не учитывать этот недостаток системы с точечными диполями, то для каждой конфигурации системы необходимо решать итеративно задачу на самосогласование поляризуемых диполей [48], что может существенно замедлять МД симуляцию (за медление расчетов МД при введении поляризуемого члена происходит всегда, но, в зависимости от метода, это замедление может быть разным). При введении подобных неаддитивных поляризуемых членов также необходимо проводить полную перепараметризацию изначальной модели, а не заимствовать параметры (например, торсионные) из старой аддитивной модели [48, 82], так как данные члены могут существенно влиять на «баланс сил», в том числе и в торсионных взаимодействиях. В недавнее время, все научные группы, которые создали популярные СП AMBER, OPLS-AA, CHARMM попытались создать поляризуемые варианты своих силовых полей. В случае СП AMBER, это простая точечно-дипольная модель, о которой кратко рассказано в начале данного раздела [48, 82]. Наследником СП OPLS-AA является поляризуемое СП PFF. Оно, как и поляризуемое СП AMBER, построено на точечных диполях без каких-либо модификаций во избегание поляризуемой катастрофы. В данном случае создатели СП полагаются на то, что обменно-отталкивающая часть ван-дер-Ваальсового взаимодействия не допустит очень близкого сближения поляризуемых диполей. Для параметризации данного СП были активно использованы квантово-механические данные, но для тестовой симуляции белка BPTI в воде [83] данное силовое поле несколько раз немного модифицировалось, чтобы улучшить описание данного белка.
При использовании специальных алгоритмов вычисления электростатической и поляризуемой части скорость вычислений падает всего лишь в 1.5 раза по сравнению с классическим СП OPLS-AA [83]. Исходное СП CHARMM было разделено на несколько ветвей в своем дальнейшем развитии. Первая представлена СП FQ-CHARMM [84, 85], в которой неаддитивная часть реализована не через точечные диполи, а через флуктуирующие заряды (FQ) с помощью фиктивной лагранжевой динамики значений отрицательных точечных «электронных» зарядов на атомных центрах (заряды положительных ядер оставались постоянными). В данной системе динамика зарядов определяется специально подогнанными параметрами электронегативности и параметрами атомной жесткости для каждого типа атома, а также лагранжевым условием на электронейтральность молекул, что приводит к эффективному перемещению зарядов по молекуле, таким образом описывая реакцию электронной плотности на внешнее поле. Этот метод удобен тем, что уравнения для зарядов выглядят так же, как и уравнения движения. При сравнении средней молекулярной поляризуемости видно, что для систем с явным квадруполем (бензол, имида-зол, формамид), который в данном СП не учитывается даже косвенно, предсказание СП недооценивает экспериментальные поляризуемости на 16%, 23.5%, 19.5%, соответственно. Новые параметры были получены с помощью КМ метода функционала плотности (DFT), который использовался для вычисления свойств молекул и димеров, что могло существенно ухудшить качество СП из-за несовершенства метода (см. раздел 1.4). Несмотря на несовершенство данной модели, она лучше описывает 6 белков, чем старая аддиттивная модель CHARMM22 [53], если судить по отклонению от нативной экспериментальной структуры при симуляции МД в течении 2 не [85]. Вторая поляризуемая ветвь [18, 86-88], созданная на основе CHARMM, базируется на так называемых осцилляторах
Друде, которые формируются из отрицательных зарядов, гармонически привязанных к положительным ядрам, таким образом формируя статический заряд в отсутствии внешнего поля, действующего на атом. При наличии внешнего электрического поля данные подвижные отрицательные заряды q отклоняются в потенциале U(x) = kx2 /2 на расстояние х от положения положительного ядра (так что атом обладает поляризуемостью a = q21 к, где q - отрицательный подвижный заряд, к - подогнанный параметр жесткости гармонического осциллятора). К сожалению, на сегодняшний день данная модель не доведена до полной биомолекулярной параметризации и не может быть применена ни к белкам, ни к нуклеиновым кислотам. Третья ветвь СП CHARMM реализована в СП Polarizable Intermolecular Potential Function (PIPF) [89] и базируется на модели Толе [81, 90]. В данной мо дели взаимодействие точечных диполей специальным образом модифицируется, чтобы избежать поляризационной катастрофы. Модификация связана с тем, что диполи представляются не в виде точечных объектов, а состоящими из зарядов, которые представляются в виде размытых экспоненциально-спадающих облаков с плотностью заряда р(и) = (a3 / S7r)e"au, где а — эффективная ширина облака, и — расстояние от ядра. В данной модели классический тензор Tlass = 3r(r,r 5 - 6/yr-3 для, создаваемого диполем ц , электрического поля Е, = Т /а"ц7 (в модели Толе Е(=Ти ) в точке г (г = г, б - символ Кронекера) заменяется на выражение Ту = Т (\-(\ + аг + а2г2/2)е аг)-а\г/-2е-аг/2 [81].
Силовое поле PIPF, как и предыдущая модель, доведено всего лишь до предсказания свойств жидкостей с алифатическими углеводородами и амидами. Однако у модели Толе есть преимущество, которое делает ее чрезвычайно привлекательной. Ее функциональная форма позволяет использовать быстрые алгоритмы вычисления электростатических взаимодействий (в том числе и поляризуемой ее части) [69, 91], подобно методу Particle Mesh Ewald (РМЕ) [33, 92, 93], который является стандартом для использования в обычных аддитивных электростатических частях СП. Создатели СП AMOEBA [69, 79, 94, 95], которое также построено на основе модели Толе [90], показали на примере многих конформаций ряда молекул одинаковой длины (от гексана с последовательным добавлением полярных элементов до цвиттер-иона аминокислоты), что СП, построенное на основе малых молекул, не хуже, а по некоторым критериям лучше, чем СП, построенное на основе различных конформаций больших, изучаемых в этой работе молекул [95]. Но главным результатом работы [95] была демонстрация факта, что неучет квадру-польных моментов на атомах равен приблизительно неучету эффекта поляризуемости и понижает средний коэффициент корреляции между предсказаниями СП и квантово-механическими данными (КМД) для всех конформаций полярных молекул с 0.9 до 0.8. При этом полный неучет всех электростатических эффектов (то есть учет только ван-дер-ваальсового отталкивания и притяжения) дает среднюю корреляцию 0.6. Из вышесказанного можно сделать вывод, что для кор
Функциональная форма СП QMPFF3
Как и в модели Друде (см. раздел 1.2.3), в QMPFF атом представляется как сумма положительно заряженного точечного ядра с отрицательным электронным зарядом, который может двигаться в гармоничном потенциале. В отличие от этой модели, в QMPFF электронная компонента моделируются не точечным зарядом, а размазанным электронным облаком, что является физически более обоснованным приближением [120]. Данное описание электронной плотности позволяет реалистично описывать так называемый эффект взаимного проникновения электронных плотностей молекул, чего не могут описать СП с точечными зарядами и мультиполями на атомных центрах [22]. В предыдущих версиях QMPFF1,2 [63, 64] электронная плотность"для каждого атома моделировались двумя электронными облаками. В отличие от этих версий QMPFF1,2 [63, 64], в QMPFF3 [22, 65] на каждый атом приходится по одному электронному облаку, которое, в общем случае, анизотропно, и изменение его положения и формы из-за изменения дипольной поляризуемой компоненты облака моделирует электронную поляризацию. Данное усовершенствование, по сравнению с ранними версиями, приводит к более экономичным вычислительных схемам в многочастичных приложениях (МД моделирование белков с QMPFF3 СП работает всего лишь в 5 раз медленнее, чем неполяризуемые СП). Таким образом, в QMPFF3 зарядовая плотность электронного облака атома а записывается в форме: Здесь Ra - позиция ядра атома, wA - его эффективная ширина, a qa — заряд, а Da(Ra)- описанный ниже (в (1.5.2) дифференциальный оператор). Заряд электронного облака записывается в форме qa = —ZA + qBA , где Z,A — be{a) положительный заряд ядра для атома а, который имеет тип А, [а] - множество атомов, химически связанных с атомом а, и qBA - параметр переноса заряда по связи В А (предполагается равенство qeA ЧАВ, ЧТО обеспечивает сохранение заряда). С помощью такого определения заряда электронного облака соблюдается электронейтральность молекулы qa + Za = 0, где MOL - множество всех aaMOL aeMOL атомов в молекуле.
Дифференциальный оператор, введенный в определении зарядовой плотности (1.5.1), действующий на сферически-симметричное облако, создает в результате анизотропное облако: ґ л V Р5 о»—, (1.5.2) 0 Ж, Вектор ta, характеризующий диполь атома а, является суммой постоянной и поляризуемой частей ta = tf + tpal. Постоянная часть описывается как tpaer - V t 5 где tab = tABnab, причем паЬ - единичный вектор, направленный Ье{а) вдоль связи от атома а к Ъ . Параметр связи tAB представляет собой постоянную «химическую» поляризуемость вдоль связи ab (стоит отметить, что не обязательно tj&fitBA)- Поляризуемая часть атомного дипольного момента tpJ = tjaxra выражается через ,ах - заданный параметр максимально-допустимого смещения облака для данного типа атома и ха -вектор смещения электронного облака, не превышающий по модулю 1, который участвует в поляризационной компоненте энергии (1.5.6). Схематическое изображение дипольной части атома а представлено нарис. 1.1. индексы а,/3 = x,y,z (Ьар- символ Кронекера, где ЬаВ = 1 при а = рк Ьав = 0 при аФ/3). Базируясь на представлении (1.5.1-1.5.2), потенциал pab ES взаимодействия между двумя зарядовыми распределениями атомов а и b отделенными расстоянием RoA=R0-R6, может быть выражен в форме (pab(Rab,wA,wB,taitb, Qa, ob) = DaDbq?(b (Rab,wA,wB) где точное выражение взаимодействия двух экспоненциально-спадающих электронных облаков с единичными зарядами и с ширинами wA, wB известно как [120] : Вычисление полной электростатической энергии между двумя анизотропными электронными облаками qaqb(pab(RabiwA,wB,ta)tb,G)a,(ob), включающей все взаимодействия вплоть до квадруполь-квадрупольных, детально выписано в [22].
Взаимодействие отрицательных электронных облаков с точечными положительными ядрами выписано в [64, 145]. Точечные ядра взаимодействуют между собой по классическому кулоновскому закону. В итоге, полная электростатическая энергия взаимодействия между двумя атомами выражается через параметры атомов: иаь =U b (Rab,Za,Zb,qa,qb,wA,wB,ta,tb,(»a,G b). В свою очередь, параметры атома а (аналогично и для атома Ь) суть следующие: qa,ta, a - определяются валентным окружением атома (qAB, tAB, QAB) И поляризуемо-деформируемым облаком, смещение которого определяется постоянным параметром іАж и подвижным вектором та, координаты которого, в свою очередь, определяются поляризационным потенциалом (см. раздел 1.5.5). Энергия обменно-отталкивающего взаимодействия UEX записывается сход ным образом с электростатическим взаимодействием (ES): u {Rab) = CACCBXXab(Rab)y гДе с" и СР сил овью параметры для атомных типов А и В, а обменный потенциал учитывает отталкивание анизотропных облаков а и Ь\ ХаЬ{ аЬ,иа ЩЛаЛь, а, »ь) = DaDaT{ аъ иа иь) гДе Два сферических диффузных облака отталкиваются по закону:
Оптимизация энергии кристалла
Тип атома», определяющий его свойства в различных взаимодействиях, зависит не только от химического сорта этого атома, но и от его ковалентного окружения. Принятая в СП ENCAD [30] и используемая нами типизация атомов для невалентных взаимодействий приведена в таблице 2.2 (при этом мы ограничились только атомами, встречающимися в белках), а более детальная типизация атомов для валентных взаимодействий может быть найдена в таблице П1.1 (Приложение 1) В качестве парциальных зарядов всех атомов всех молекул были использованы принятые в силовых полях, рассчитанные на основе квантовой механики, «RESP заряды» [49, 112, 143], дающие оптимальную картину пространственного распределения электростатического потенциала молекулы. Квантовую оптимизацию конфигураций молекул проводили, стартуя с их кристаллической структуры, на основе рекомендуемого [49] набора базисных функций HF/6-31G (см. [121]) с помощью программы PC GAMESS/Firefly QC [121], которая частично базируется на коде GAMES S (US) [170]. Получив оптимальную электронную конфигурацию молекулы, заряды ее электронных облаков приписывались ядрам ее атомов утилитой resp (из пакета Antechamber 1.27 [171]). Для получения параметров невалентного взаимодействия необходимо оптимизировать кристаллические ячейки и отдельные молекулы в газе (см. раздел 2.3). Аналитически вычисляемый градиент энергии по координатам атома и параметрам ячейки позволяет использовать быстрый алгоритм минимизации LBFGS ("Limited memory Broyden-Fletcher-Goldfarb-Shanno method") [172]. Энергия молекулы в кристалле: [Um «ооо + noncovWooo] + Уг КЕ" (Мооої ю,а + т2Ъ + т3с)] получается делением суммы всех взаимодействий внутри одной его (кристалла к) кристаллической ячейки, координаты атомов в которой задаются векторами {r}000, и полусуммы всех невалентных взаимодействий атомов этой ячейки с лежащими в пределах Яс от них (см. уравнение (2.1.4)) атомами других ячеек (сдвинутыми от {г}000 соответственно векторам а,Ь,с решетки кристалла к) на число Z молекул в ячейке кристалла. В предыдущих работах по получению параметров СП из кристаллов [41, 56, 57, 60] оптимизация структур кристаллов для каждого набора параметров невалентных взаимодействий либо не проводилась вообще, либо относилась только к параметрам кристаллической решетки, а положения молекул в рамках кристаллической ячейки были заморожены. Рост вычислительных мощностей к настоящему времени позволил нам впервые провести полную оптимизацию структур кристаллов.
Однако у кристаллической ячейки обычно есть внутренние симметрии, которые, в принципе, не должны нарушаться в ходе оптимизации. Если учесть эти симметрии, пространство оптимизируемых координат сужается, и пропорционально, в алгоритме минимизации падает число итераций (хотя количество вычислений парных взаимодействий в одной итерации не уменьшается). В базе кристаллов CSD [160] даны: векторы acr, bcr, ссг ячейки кристаллической решетки; заданные в базисе {ac/.,bcr,cc/.} координаты х пе\...,х1ег всех А атомов молекул в ассиметричной единице ячейки и данные по внутренней симметрии ячейки. Последние заданы в виде 3x3 матриц ортогональных преобразований О c=I s (которые и в косоугольных ячейках не нарушают структуру молекул) и соответствующих векторов смещения rij[[ s (где S - число используемых симметрии). Минимизация энергии ведется в пространстве внутренних координат ячейки. На вход подаются координаты г атома к, соответствующие 1-й симметрии ячейки: r ke =(x;y\z) x; этой 1-й симметрии отвечает 0 =1- единичная матрица преобразований координат, и rff = (0;0;0). Затем по симметриям =2, ..., S вычисляются внутренние координаты атома к во всех остальных молекулах ячейки: г ктГ = О,r k mer + rfift (где О - уже не единичные матрицы); после чего можно получить декартовы координаты атомов к во всех молекулах ячейки 0 Ь2 с2 а преобразования А = По декартовым координатам считается энергия U и силы Carl dU дг После этого силы, с учетом симметрии, конвертируются в силы во внутренних координатах формулы энергии кристаллической ячейки или молекулы. Часть сил, связанных с валентными взаимодействиями и действующих на атом /, является производной функ цией валентных взаимодействий по вектору rfart и является функцией координат атома и его валентных (вплоть до торсионных ±3 атома по цепи) соседей dUmv/dr =Fmv(r \rfj\.r ). Детальное описание производных валентной части взаимодействия находится в файле http://phys.protres.ru/resources/FFS/Al.doc. Невалентная часть сил, действующая на атом / в центральной ячейке (0,0,0), складывается из суммы действия атомов j внутри ячейки, что соответствует члену noncovWooo в формуле (2.2.3) Для взаимодействия между ячейками У2 2L t C ((r)ooo5wia + Щр + пус)] т1,т2,гщ эффективная сила, действующая на атом і, состоит из действия атомау в ячейке п и из такого же действия на образ атома і в противоположной ячейке (относительно ячейки (0,0,0)) образом атома у в центральной ячейке, поэтому результирующая сила, действующая на атом /, должна записываться без коэффициента И ,
Параметры вспомогательного силового поля FFSubl и тестирование их качества
Расчет «вакуумных» потенциалов и их тестирование является отправной точкой работы, позволяющей проверить метод получения параметров. Потенциалы этого «вакуумного» СП (поскольку оно основано на энергии сублимации, назовем его FFSubl) рассчитывались двумя методами: с соблюдением кристаллографических симметрии в ходе локальной минимизации энергий структур молекул в кристаллах (см. выше) и без такого учета. Таблица 2.3 показывает, что полученные в обоих случаях параметры СП примерно одинаковы; однако общая функция невязки А (между вычисленными и экспериментально определенными свойствам кристаллов, см. формулу (2.3.5) и предыдущие формулы (2.3.1) - (2.3.4)) существенно ниже в случае минимизации с соблюдением кристаллографической симметрии. При этом основной выигрыш достигается за счет лучшего описания геометрии ячейки кристалла. По-видимому, это объяснятся тем, что оптимизация с учетом симметрии работает в пространстве с меньшим числом параметров и ее более «усредненный» путь не застревает в малых «случайных» локальных минимумах энергетической поверхности. Предсказанные СП FFSubl величины A/cryst-vap как с учетом симметрии, так и без него, достаточно хорошо согласуются с экспериментом: коэффициент корреляции г = 0.954 в первом случае и г=0.947 во втором. Частные функции невязки А(=123 для каждого кристалла оптимизированного с учетом симметрии представлены в таблице П1.6 (Приложение 1). Полезно также сравнить результаты работы лучшего (полученного с учетом симметрии кристаллов) из наших СП FFSubl с результатами работы других СП, для которых в научной литературе можно найти соответствующие данные. Сравнение (таблица 2.4) проводилось с СП CFF II [41] (подогнанным под кристаллические структуры без учета симметрии в них) и с СП QMPFF3 (подогнанным под квантово-механические расчеты взаимодействий простых молекул, и также не учитывающим симметрии кристаллов) [65]. В таблице 2.4 сравниваются относительные точности воспроизведения энергий взаимодействий и объемов кристаллических ячеек по пересекающимся (для данной работы и для работ [41, 65], соответственно) множествам кристаллов. Из таблицы 2.4 видно, что СП FFSubl (построенное как с учетом, так и без учета симметрии ячеек кристаллов) лучше воспроизводит и энергии молекул в кристаллах, и объем их ячеек. Это может быть связано с более тщательной процедурой оптимизации, применяемой в настоящей работе.
Однако этот эффект может быть отчасти связан и с тем, что при оптимизации CFF II и QMPFF3 для экспериментальной оценки энергии выхода молекулы из кристалла (АС/ ) в ра ботах [41, 65] было использовано чуть более грубое, чем при построении базы данных для FFSubl, приближение ЛС/ехр = -АЯ5иЫ -2RT (сравни с формулой (2.1.7) и примечаниями к ней). И еще этот эффект может быть связан с тем, что параметры СП в FFSubl подгонялись под одну температуру 25 С, к которой экстраполировались все экспериментальные данные, тогда как в СП CFF II параметры подгонялись под величины ДС/ехр, частично экстраполированные к 0 К, а частично нет, а в СП QMPFF3 воспроизводилась сублимация, происходящая при разных температурах. 2.4.2 Параметры «водного» силового поля FFSol и тестирование их качества Получение «водных» потенциалов СП FFSol является одной из целей данной диссертационной работы. В таблице 2.5 даны атомные параметры невалентных взаимодействий для «водного» силового поля FFSol, соответствующие минимизации общей функции невязки А на полном наборе 58 используемых молекул при растворении их кристаллов в воде. Для контроля там же даны аналогичные параметры, полученные при минимизации А на двух половинах этого набора. Наборы выделялись таким образом, чтобы каждый из них представлял приблизительно половину молекул для атомов каждого типа. Первый из наборов состоит из 28 органических молекул, а также молекулы воды, а второй — из остав шихся 29 органических молекул и молекулы воды (которая входит в оба набора, так как только она включает атом типа W). Во всех случаях минимизация проводилась с учетом симметрии кристаллов, в соответствии с опытом, полученным нами при получении «вакуумных» параметров силового поля FFSubl (таблица 2.3). Частные функции невязки A,=12i3 вычисленные с помощью параметров полученных на полном наборе для каждого кристалла оптимизированного с учетом симметрии представлены в таблице Ш .6 (Приложение 1). Сравнение таблиц 3.3 и 3.4 показывает, что параметры «водного» СП FFSol довольно сильно отличаются от параметров «вакуумного» СП FFSubl. При этом основная разница связана с ожидаемым резким (по сравнению с FFSubl) ростом є (и, соответственно, ослаблением электростатических взаимодействий) в «водном» СП FFSol и с уменьшением в нем параметров В , отвечающих притяжению полярных, особенно водородных (D) атомов, при некотором росте параметров 5 для неполярных атомов. Интересно, что в "вакуумном" СП FFSubl оптимизированная величина є лежит между типичной величиной диэлектрической проницаемости кристаллов органических молекул ( 3, см. [115, 162]) и диэлектрической проницаемостью вакуума (1.0), а в "водном" СП FFSol є лежит между диэлектрической проницаемостью кристаллов органических молекул и воды («80). Таблица 2.5 показывает также, что получаемые параметры СП FFSol обычно мало зависят от того, на каком наборе проводилась их оптимизация. Это свидетельствует о «переносимости» СП FFSol с одних молекул на другие. Кроме того, таблица 2.5 показывает, что величина функции невязки А, полученная на всем наборе молекул, мало отличается от полученной на двух «половинных» наборах.
При этом сравнение таблиц 2.3 и 2.5 показывает, что функция невязки для «водного» случая вдвое больше, чем для вакуумного, что отвечает увеличению погрешностей приблизительно на 40% (см. величины АЕ/Е в таблице 2.6). Рост невязки в СП FFSol естественен, так как наша модель, неявно описывающая водное окружение, заметно менее строга, чем модель, описывающая взаимодействие в вакууме. Гораздо более любопытно то, что этот рост относительно невелик, — это показывает, что используемая нами модель влияния водного окружения далеко не плоха. Замечательно также, что оба наших сравнительно простых поля, FFSol и FFSubl, воспроизводят экспериментально наблюдаемые энергии и объемы кристаллов ничуть не хуже, чем значительно более детальные силовые поля CFF II и QMPFF3. Таблица 2.6 показывает, что и у использующихся «вакуумных» СП CFF И, QMPFF3, нового «вакуумного» СП FFSubl, и у нового «водного» СП FFSol относительная погрешность расчета энергии составляет 6.3 - 9.7%, а относительная погрешность расчета объема кристаллической ячейки — 2.8 — 4.1%. Рисунок 2.1 показывает, как силовые поля FFSol (в том числе полученные из «половинных» наборов молекул) описывают перенормированную водой разность «потенциальных энергий» основных состояний молекулы в кристалле и в воде, AC/cryst.aq. На этом рисунке для СП FFSol можно увидеть, что лишь молекула № 49 (фенотиазин) находится далеко от диагональной линии. Большое отклоне
