Содержание к диссертации
Введение
II. Литературный обзор 5
II. 1. Заряды в молекулярном моделировании 5
II.2. МЭП как фундаментальное свойство 7
ІІ.3. Типы зарядовых схем 11
II.3.1. Неэмпирические методы расчета электронной плотности 11
II.3.2. Полуэмпирические методы расчета электронной плотности 15
II.3.3. Эмпирические методы 16
II.3.4. Заряды, воспроизводящие различные свойства 30
II.4. Топологическая симметрия зарядов 36
II.5. Сравнение существующих схем 37
II.6. Выводы из анализа литературы 39
III. Теоретическая часть 41
III. 1. Основания и постановка задачи 41
ІІІ.2. Динамическая релаксация ЭО 41
III.3. Полное выравнивание ЭО с топологически симметричной функцией энергии 50
III.4. Сравнение с существующими методами 52
ІІІ.5. Практические схемы 56
III.5.1. Модель молекулярного графа 57
III.5.2. Модель орбитального графа 57
III.5.3. Модель динамической релаксации ЭО 58
III.5.4. Модель полного выравнивания ЭО с топологически симметричной функцией энергии системы 59
IV. Детали вычислительного эксперимента 60
IV. 1. Постановка задачи 60
IV.2. Метод 60
IV.2.1. Выборка структур 60
IV.2.2. Оптимизация геометрии 61
IV.2.3. Расчет МЭП на решетке точек 61
IV.2.4. Целевые функции 62
IV.2.5. Использование атомных типов 64
IV.2.6. Методика оптимизации 67
IV.2.7. Контроль качества 70
IV.2.8. Программный комплекс 70
IV.2.9. Расчет зарядов, используемых для сравнения 71
V. Обсуждение результатов 73
V. 1. Воспроизведение RESP-зарядов зарядами схем МГ и ОГ 73
V.2. Воспроизведение МЭП зарядами схем МГ, ОГ, ТСФЭ и ДРЭО 83
V.2.I. Общие особенности этапа 83
V.2.2. Схема МГ 87
V.2.3. Схема ОГ 92
V.2.4. Схема ТСФЭ 96
V.2.5. Схема ДРЭО 101
V.2.6. Выводы этапа оптимизации 106
V.3. Общий анализ результатов 106
V.3.I. Сравнение эмпирических схем расчета зарядов 106
V.3.2. Особенности воспроизведения МЭП 112
V.3.3. Сравнение с литературными схемами 117
V.3.4. Анализ сложностей воспроизведения МЭП атомными зарядами 133
V.3.5. Примеры расчета зарядов и МЭП 142
V.3.6. Дальнейшее развитие 144
VI. Выводы 150
VII. Список литературы 151
Приложение 164
- Неэмпирические методы расчета электронной плотности
- Полное выравнивание ЭО с топологически симметричной функцией энергии
- Использование атомных типов
- Анализ сложностей воспроизведения МЭП атомными зарядами
Введение к работе
Межмолекулярные взаимодействия играют важнейшую роль как при исследовании живых систем на молекулярном уровне, так и при разработке новых перспективных материалов с заданными свойствами. Такие исследования - через направленный органический синтез - приводят к соединениям, имеющим широкое практическое применение и большое значение в фармацевтической химии, медицине, промышленности и сельском хозяйстве. Атомные заряды оценивают наиболее важный — электростатический — вклад в межмолекулярные взаимодействия, поэтому они нашли широкое применение в различных областях молекулярного моделирования. Исторически атомные заряды применялись для предсказания относительной химической активности органических соединений. К сожалению, атомные заряды невозможно определить уникальным образом, в связи с чем возникло большое разнообразие методов их расчета. Опыт применения атомных зарядов в различных областях молекулярного моделирования (молекулярная динамика, докинг, конформационный поиск, корреляции «структура-свойство»/«струк-тура-активность» QSPR/QSAR) сформировал целый ряд в чем-то противоречивых требований к зарядам и схемам их расчета: высокая скорость расчета, хорошее воспроизведение электростатического потенциала вокруг молекул (молекулярного электростатического потенциала - МЭП), соблюдение топологической симметрии зарядов для интенсивного исследования конформационного пространства. Предложенные к настоящему моменту схемы не удовлетворяют части этих требований. Кроме того, во многих случаях теория, положенная в их основу не вполне прозрачна физически, что препятствует их последовательной модификации для достижения оптимального соотношения точности описания и скорости вычисления зарядов - важного дополнительного требования при разработке перспективных методов многоуровневого моделирования, действующих на масштабах от небольших органических структур до наноразмерных ассоциа-тов макромолекул.
Таким образом, существующая потребность в эмпирических схемах расчёта атомных зарядов, удовлетворяющих всем вышеописанным требованиям, делает несомненной актуальность задач, связанных с разработкой подобных методов расчёта.
Неэмпирические методы расчета электронной плотности
Наиболее известным методом этой группы является анализ орбитальной заселенности по Малликену [19], в котором заряд на каждом атоме складывается из значения заряда ядра и суммы заселенностей орбиталей, центрированных на данном атоме. Краеугольным камнем этого метода является способ разделения вкладов существенно перекрывающихся орбиталей, участвующих в описании валентных электронов. В данном подходе вклад от перекрывающихся орбиталей, центрированных на различных атомах, делится поровну. Несмотря на значительную критику, метод анализа заселенности по Малликену реализован в каждом квантово-химическом пакете, поскольку не требует значительных вычислительных затрат после того, как произведено самосогласование электронной плотности, и не использует каких-либо дополнительных допущений. Вариантом анализа заселенности является также схема, предложенная Лёвдиным [20].
Изначально оба метода применялись для анализа результатов расчета, основанных на самосогласованном решении уравнений Хартри-Фока. При этом была продемонстрирована существенная зависимость величин зарядов от используемого базисного набора. Величины дипольного и квадрупольного моментов структуры, оцененные при помощи этих зарядов, часто заметно отличаются от величин, рассчитанных непосредственно из волновой функции [109].
Анализ натуральных орбиталей связей (NBO) [21, 22] позволяет выделить компактную и привычную для химика картину образования связей и взаимодействий на молекулярном уровне, особенно полезную для анализа донорно-акцепторых ван-дер-вааль-совых комплексов [23] и предлагает свое определение атомных зарядов в рамках натурального анализа заселенностей (NPA) [22]. Учет электронной корреляции при получения атомных зарядов, по всей видимости, не оправдан из-за неоправданно высоких вычислительных затрат, поскольку при переходе от распределенного пространственного описания зарядовой плотности к ее дискретному представлению (набор зарядов) происходит колоссальная редукция информации. Возможно в некоторых специальных случаях оправдано использование теории возмущения Меллера-Плессе второго порядка (МР2), ввиду относительной важ ности вносимых поправок и распространенности эффективных реализаций в квантово-химических пакетах. Альтернативным подходом к расчету зарядового распределения является подход теории функционала электронной плотности (ТФП, DFT) [24], основанный на решении уравнений Кона-Шема [25, 26]. Согласно теореме Хохенберга-Кона [27, 28] распределение электронной плотности в пространстве вокруг ядер системы в основном состоянии минимизирует полную энергию этой системы, определенную как функционал электронной плотности. Подход ТФП интересен тем, что теоретически может привести к абсолютно верному решению при условии, что точно задан вид так называемого обменно-корреляционного функционала [24], который к настоящему моменту не известен, хотя предложено несколько его аппроксимаций [29, 30], получивших широкое распространение, а известные на данный момент функционалы уже достаточно точны для многих химических приложений [31]. В рамках ТФП электронная корреляция учитывается до некоторой степени даже в самом приближенном вычислении, однако мера учета остается неизвестной. Недостатком ТФП является то, что оценка энергии системы не является оценкой сверху, как это выполняется для традиционных методов квантовой химии, т.е. не соблюдается вариационный принцип. Достоинством ТФП в приближенных вычислениях в достаточно больших системах является лучшее масштабирование по сравнению с квантово-химическими методами, особенно с использованием локальных обменно-корреляционных функционалов [31].
Вычисление зарядов из результирующей электронной плотности системы осуществляется аналогично описанному выше способу. Другой способ получения значений атомных зарядов исходя из результатов квантово-химических расчетов состоит в разбиении электронной плотности молекулы на атомные составляющие. Этот подход основан на свойстве атома во многом сохранять свою идентичность при образовании молекулы. Идеи, представленные в пионерской работе Политцера и Харриса [32], в которой рассматриваются простые линейные молекулы, оказалась плодотворной и была развита в последующих работах в двух направлениях: методе Хиршфельда и теории атомов в молекуле Бейдера. Элегантная схема разделения («паевых вкладчиков») электронной плотности была предложена Хиршфельдом [33]. В данном подходе электронная плотность молекулы разделяется на атомные вклады пропорционально весам, которые имеют сферически усредненные электронные плотности данных атомов в гипотетической «промолекуле», образованной невозмущенными атомными плотностями. Электронная плотность промолекулы и весовой вклад («пай») соответствующего атома представляются в виде где pi"(r) - электронная плотность невозмущенного атома / и w,(r) - вклад электронной плотности атома / в электронную плотность промолекулы в точке пространства г. Когда весовые вклады w,(r) определены, электронная плотность, принадлежащая данному атому, принимается пропорциональной доле электронной плотности молекулярной системы pmolir) (рассчитанной или полученной из экспериментальных данных) где рї"/то1(г) - электронная плотность, приписываемая атому / в данной молекуле. Частичный заряд на атоме получается интегрированием электронной плотности-, принадлежащей данному атому по всему пространству где q{ - частичный заряд на атоме /, Z, - заряд ядра атома /. Позднее было показано, что подход эквивалентен минимизации функционала дефекта энтропии (минимальной потери информации) [34]. Дальнейшее развитие метода привело к появлению «итеративного метода Хиршфельда» [35]. Как показывают авторы [36], этот подход соответствует минимизации того же функционала, с ограничениями, отличными от ограничений «метода пайщиков». Итеративная схема исправляет основной недостаток оригинальной схемы, состоящий в зависимости зарядов от способа образования промолекулы из отдельных «атомных» компонентов.
Полное выравнивание ЭО с топологически симметричной функцией энергии
Экспериментальные данные и химический опыт подтверждают существенное сохранение свойств (в том числе зарядового распределения) фрагментов молекул в различных химических окружениях, что положило основу классификации соединений в органической химии на основании наличия определенной функциональной группы. На основании этой идеи предложен аддитивный подход к расчету частичных атомных зарядов из значений так называемых связевых зарядовых инкрементов, определяющих величину переноса заряда, осуществляемого вдоль связи от одного атома к другому при гипотетическом образовании молекулы. Конечный заряд на атоме получается суммированием (с учетом знака) величин зарядовых инкрементов по всем связям, в которых участвует данный атом. Фактически такой подход эквивалентен одному шагу итераций метода Гастайгера. Такая модель использована при разработке зарядовых схем силового поля MMFF94 [89, 90, 91], в котором зарядовые инкременты были оптимизированы для воспроизведения дипольных моментов молекул, а также энергий взаимодействий и равновесных геометрий комплексов с водородной связью.
Основными достоинствами зарядовых инкрементов является скорость расчета зарядов, автоматическое соблюдение топологической симметрии и значительная гибкость в представлении зарядовых распределений самых разнообразных органических соединений, что важно в исследованиях, связанных с фармацевтической и медицинской химией. Обратная сторона гибкости - необходимость значительного числа параметров (зарядовых инкрементов), растущего квадратично от числа используемых атомных типов, поскольку каждый зарядовый инкремент зависит от пары атомных типов. Зарядовые инкременты автоматически сохраняют полный заряд молекулы, поэтому структуры, несущие формально заряженные группы также могут быть представлены в рамках данного подхода. В этом случае начальное значение заряда на некоторых атомах будет ненулевым. Заряды на основе зарядовых инкрементов (без введения специальных атомных типов) не могут описать индуктивный эффект, поскольку по определению отражают перенос заряда только между непосредственно связанными атомами. По этим же причинам, возникают сложности в структурах, где эффекты делокализации электронной плотности (неаддитивные эффекты) значительны, например, в сопряженных и ароматических системах. Несмотря на указанные недостатки, метод зарядовых инкрементов позволяет получать заряды, воспроизводящие электростатические свойства с точностью, достаточной для применения в классических силовых полях молекулярной механики (каковым является MMFF94).
В качестве развития подхода была проделана работа по параметризации зарядовых инкрементов для воспроизведения квантово-химического X /6-31G(d) электростатического потенциала одновременно вокруг выборки разнородных структур [92]. Выявленные закономерности, а также принятая двухстадийная методика использования зарядовых инкрементов из работ Крамера (см. П.3.4.6) заложили основу для дальнейшей работы по созданию зарядовой схемы АМ1-ВСС [93], в которой заряды, получаемые полуэмпирическим методом AMI, корректируются зарядовыми инкрементами для лучшего воспроизведения расчетного МЭП молекул. При помощи данного подхода не только удалось достичь качества описания МЭП, лишь немногим уступающего зарядам RESP [94, 95] (непосредственно воспроизводящим электростатический потенциал для каждой структуры по отдельности), но и добиться большей устойчивости зарядовых схем, что необходимо для корректной работы силовых полей, для которых эти зарядовые схемы разрабатывались. Использование в качестве стартовых AMI-зарядов позволяет решить указанные выше проблемы с учетом индуктивного и мезомерного эффектов, а также распределения заряда в формально заряженных молекулах и структурах с существенной делокализацией электронной плотности. Значения зарядов при этом теряют топологическую симметрию, т.к. отвечают одному из конформеров (в том числе вращательных), найденному в результате оптимизации методом AMI. Кроме того, полуэмпирический расчет не всегда выполняется достаточно быстро по меркам требований высокой производительности, предъявляемых в области медицинской химии - основного потребителя результатов молекулярного моделирования.
Распространенной практикой в моделировании при помощи силовых полей биополимеров, как соединений, составленных из небольшого числа различных мономеров (аминокислот - в белках, нуклеотидов - в нуклеиновых кислотах), является использование табулированных значений зарядов составляющих мономеров [17]. При этом предполагается, что перенос заряда от мономера к мономеру отсутствует. В противном случае табличные значения не суммировались бы к целым числам, отвечающих зарядовым состояниям биомолекул. Используя стандартные значения зарядов, разработчики силовых полей получают возможность добиться их лучшего согласования с ван-дер-ваальсовыми параметрами для совместного использования. Такой подход, однако, не применим для произвольных органических соединений, поскольку, во-первых, требование нулевого переноса заряда между табулированными фрагментами оказывается слишком строгим, и, во-вторых, число фрагментов, необходимое для описания произвольной органической структуры должно быть достаточно большим.
Для получения топологически симметричных и одновременно зависящих от химического окружения значений зарядов был предложен оригинальный подход [96, 97, 98], в котором молекулярная структура сравнивается с замкнутой электрической цепью, ЭО атомов - с источниками электродвижущей силы (э.д.с), химические связи между атомами - с электрическими сопротивлениями, зависящими от кратности связи. После замыкания цепи устанавливаются стационарные токи и соответствующие им падения напряжения на дополнительных сопротивлениях (рис. 2.1), которые соответствуют установившимся в системе эффективным значениям ЭО.
Использование атомных типов
Проведенный литературный обзор основных способов расчета значений зарядов с учетом требований для типовых приложений молекулярного, и в частности биомолекулярного, моделирования, дает основание для классификации зарядовых схем по набору .j критериев, выдвинутых в начале раздела П.З. В результате можно сделать ряд выводов.
Заряды, основанные на квантово-химических расчетах непрактичны в силу невыгодного соотношения большой вычислительной трудоемкости и умеренного качества воспроизведения МЭП для схем, основанных на анализе электронной плотности: как к анализу заселенностей Малликена, Левдина и др., так и к анализу топологии электронной плотности Бейдера. При вычислении значений атомных зарядов происходит значительная редукция информации, что приводит к некоторому произволу.
Для совместного использования с классическими силовыми полями молекулярной механики с точки зрения адекватности и точности лучшими являются заряды, воспроизводящие [XO/6-31G(d)] МЭП на решетке вокруг структур. Такие заряды неявным образом отражают поляризацию электронной плотности, вызванную помещением органической молекулы в водную среду. Однако эти методы не пригодны для массового расчета зарядов для больших выборок структур и отдельных молекул больших размеров, т.к. требуют качественной начальной геометрии и достаточно трудоемки. Топологическая симметрия накладывается ограничениями в процедуре МНК подбора зарядов, причем для составления ограничений часто требуется квалифицированное вмешательство пользователя. Статистическая определенность зарядов атомных центров в процедуре МНК часто сильно отличается. Заряды на слабоопределенных центрах могут принимать неадекватные величины, ограничивающие согласование зарядов с остальными параметрами силовых полей [93]. Тем не менее, именно с зарядами этого типа параметризуются классические силовые поля.
Полуэмпирические методы значительно менее требовательны к вычислительным ресурсам, чем неэмпирические, однако сохраняют основные недостатки последних. Из-за большей вычислительной доступности успешно применяются комбинированные схемы, в которых заряды, полученные в результате полуэмпирического расчета, корректируются набором эмпирически оптимизированных инкрементов для воспроизведения получаемыми зарядами заданных свойств (МЭП - для схемы АМ1-ВСС, дипольных моментов - для методов СМХ, Х=1,2,3). Заряды зависят от химического окружения, но не отражают топологической симметрии.
Эмпирические схемы привлекательны по вычислительной трудоемкости и гибкости в подборе параметров, однако для большинства из них (кроме, например, схемы MMFF94) в настоящий момент не представлена параметризация, направленная на воспроизведение МЭП. Методы, основанные на принципе полного выравнивания ЭО, требуют задания начальной геометрии, а учет влияния изменения геометрии на значения зарядов либо не всегда адекватен [90], либо трудоемок. К тому же способ учета этого влияния должен быть согласован с энергетической функцией (ММ энергией, оценочной функцией) для успешного применения. Топологическая симметрия нарушается тем же средством, которое используется для учета влияния химического окружения, - функцией энергии, зависящей от межатомных расстояний. 5. Эмпирические схемы, использующие неполное выравнивание ЭО, такие как схема Гастайгера или зарядовые схемы МГ и ОГ, по многим критериям оказываются предпочтительнее, поскольку не требуют начальной геометрии, а получаемые заряды зависят от химического окружения, обладают топологической симметрией и хорошо коррелируют со многими наблюдаемыми свойствами [73], что свидетельствует об адекватности используемого приближения. К наиболее существенным недостаткам данных схем в области количественного моделирования можно отнести параметризацию, основанную на атомарных и молекулярных электронных свойствах, и, как следствие, недостаточное качество воспроизведения МЭП. Заряды схемы Гастайгера существенно меньше по абсолютной величине [90] зарядов, адекватно воспроизводящих МЭП (Мерца-Коллмана, RESP, CHELPG и MMFF94). 6. Зарядовые инкременты схемы MMFF94 позволяют добиться адекватного качества описания МЭП [92] и большой гибкости в рамках аддитивной схемы основных вкладов в электростатический потенциал разнообразных органических структур. Обратная сторона гибкости - квадратичный рост числа зарядовых инкрементов. Недостаток аддитивной схемы - неучет индуктивного эффекта и эффекта насыщения, наблюдаемого, например, по снижению абсолютного заряда на атоме галогена при последовательном замещении водорода в метане. 7. Заряды, подобранные для расчета определенного свойства, такого как молекулярные мультиполи, энергии сольватации, ИК-спектров (GAPT) и др. часто неявным образом отражают некоторые свойства, присущие их применению и не отражают особенности воспроизведения МЭП. Кроме того, такие заряды, как правило, непригодны в потоковом высокопроизводительном моделировании. 8. Заряды, подобранные как параметры силового поля для описания эталонных энергий взаимодействия, а также среднестатистических свойств ансамблей, хотя и являются наилучшими по достижимой точности аппроксимации, требуют наибольшего привлечения вычислительных ресурсов, а также искусства экспериментатора. В силу того, что заряды в данном случае могут вместить максимально возможное количество неявной информации, и того, что заряд на атоме не является строго локальным свойством, такие заряды нередко плохо переносимы из одного силового поля в другое.
Анализ сложностей воспроизведения МЭП атомными зарядами
Рассмотрим свойства полученного нестационарного решения. Вклад компонент решения, отличных от первого собственного вектора, экспоненциально убывает со временем эффективной эволюции системы. В пределе при t—»со остается только стационарное решение. Однако поскольку собственные вектора матрицы Лапласа L зависят от молекулярной топологии (связности) системы, то и решение относительно ЭО (3.21) и зарядов (3.20) также зависит от молекулярной топологии в любой произвольный момент времени. Более того, вклады собственных векторов с большими собственными значениями, отвечающие существенным разностям ЭО на смежных вершинах-атомах графа-молекулы, убывают быстрее (с показателем экспоненты, пропорциональным собственному значению). Эти же компоненты соответствуют векторам, изменение ЭО вдоль которых приводит к наибольшему повышению энергии системы.
Итак, начав процесс с произвольного вектора ЭО, отвечающего нормировке на полный заряд системы и с начальными зарядами, не нарушающими топологической симметрии, система эволюционирует согласно уравнениям (3.21) так, что наиболее энергетически невыгодные компоненты убывают значительно быстрее. В целом система стремится к стационарному решению (3.16), отвечающему полному выравниванию ЭО, значения зарядов в котором не зависят от химического окружения. При этом процесс эволюции системы уместно назвать динамической релаксацией ЭО, поскольку значения ЭО для непосредственно связанных атомов сближаются гораздо быстрее, чем ЭО более удаленных атомов. Безусловно, функция (3.19) является очень грубым способом описания зависящей от заряда энергии системы, поэтому поиск зарядов, точно минимизирующих энергию (3.19) нецелесообразен. Вместо этого можно, начав с некоторого разумного приближения для распределения зарядов (например, #/=0 для нейтральных молекул), провести релаксацию системы, например, до тех пор, пока уменьшение энергии на шаге не достигнет ошибки в приближении энергии (3.19). Приведенные рассуждения составляют предлагаемый принцип динамической релаксации ЭО.
В более общем случае, когда значения жесткости атомов (или орбиталей) различаются (?7, ]: 1), невозможно получить удобного аналитического решения. Более того, решение полной задачи на собственные вектора налагает дополнительные требования на вычислительные ресурсы и память при проведении прикладного моделирования с использованием рассчитываемых зарядов. Более практичный подход, учитывая приближенную природу описания функции энергии системы, состоял бы в релаксации значений ЭО (и энергии) до выполнения заданного критерия. Этот процесс можно провести путем конечно-разностного интегрирования уравнений (3.12). Заменяя производные конечными приращениями в левой части уравнений (3.12), можно получить две простейшие схемы конечно-разностного интегрирования обыкновенных дифференциальных уравнений: явную (3.23) и неявную (3.24) схемы Эйлера, отличающиеся тем, с какого шага к используется заряды q(A) в правой части (3.12) где At - шаг дискретизации по времени; rjo=diag(ri\Q,..., rjNat) - диагональная матрица жесткости; /0 - вектор атомных ЭО; qw - заряды на итерации /, причем q(0 = qo.
Явный метод Эйлера вычислительно более прост - не требуется обращение матрицы или решения системы нелинейных уравнений в более общем случае. Однако в случае реальных химических структур получаются так называемые «жесткие системы» уравнений [133], требующие для сохранения точности и устойчивости применения неявных методов. В данном случае жесткость системы определяется отношением наибольшего к наименьшему собственных значений матрицы Лапласа и растет с увеличением длины ациклических частей молекулярного графа. Сохранение устойчивости в явной схеме требует применения малого шага по времени, приводящего к неоправданно детализированному (в силу грубости функции энергии) описанию динамики. Более сложные методы интегрирования нецелесообразны из-за приближенной природы самого вычисления, некоторого произвола в точке остановки релаксации, а также с целью сохранения разумного соотношения вычислительных затрат и качества получаемых результатов.
Нарис. 3.3 представлены результаты интегрирования уравнений различными схемами для молекулы ацетамида с параметрами [эВ]: хс =5.343, XN =6.899, хо =8.741, % н =4.528, TJC =10.126, лы =11.760, г0 =13.364, п0н =13.890. Нарис. 3.3а шаг интегрирования выбран относительно небольшим. Неявная (синяя) и явная (красная) схемы ведут себя аналогично точному решению (левая часть) - энергия стремится к своему предельному значению. Модуль вектора зарядов (правая часть) для явной схемы после небольшого всплеска также следует поведению точного решения. При увеличении шага (рис. З.Зб) кривая энергии для явной схемы демонстрирует заметное отклонение от точного решения, а модуль вектора зарядов - сильные осцилляции, сходящиеся в пределе для данной длины шага и расходящиеся при его дальнейшем увеличении, что приводит также к катастрофической потере точности.
