1. На главную
  2. Электронные версии
  3. Журнал физической химии
  4. T. 93, № 10, 2019

Журнал физической химии, 2019, T. 93, № 10, стр. 1477-1493

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

А. А. Рыбаков a, И. А. Брюханов b, c, А. В. Ларин a, *

a Московский государственный университет имени М.В. Ломоносова, Химический факультет
Москва, Россия

b Научно-исследовательский институт механики МГУ им. М.В. Ломоносова
119192 Москва, Россия

c Российская академия наук, Институт машиноведения им. А.А. Благонравова
101000 Москва, Россия

* E-mail: nasgo@yandex.ru

Поступила в редакцию 15.03.2019
После доработки 15.03.2019
Принята к публикации 09.04.2019

Полный текст (PDF)

Аннотация

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

Ключевые слова: методы расчета атомных мультипольных моментов, метод экстраполяции атомных мультипольных моментов в кремниевых цеолитах, электростатический потенциал, электростатическое поле

ВВЕДЕНИЕ

Повышение скорости квантово-химических расчетов и достижение линейной масштабируемости времени расчетов относительно увеличения размеров системы является одной из наиболее важных проблем вычислительной химии. Для ее решения методы приближенного расчета двухэлектронных интегралов по орбиталям удаленных друг от друга атомов часто обращаются к представлению соответствующих электронных распределений в виде орбитальных мультипольных моментов (ММ), оболочечных (внутренних и внешних) ММ и атомных ММ (АММ) [1]. В пакете GAUSSIAN, начиная с версии 16 [2], подход fast multipole method (FMM), оперирующий с другим способом деления пространства на ячейки и расчетом ячеечных ММ [3], по умолчанию используется во всех типах расчетов для уменьшения числа двухэлектронных интегралов, т.е., для той же цели, что и метод Саундерса и др. [1] в пакете CRYSTAL [4] версии 95 и более поздних. Континуальные модели “раствор–растворенная молекула” оперируют с мультипольными разложениями для описания энергии взаимодействия [5]. Как было показано десятки лет назад [6], расчет электростатической энергии через распределенные (по связям и атомам) ММ приводит к расширению точности ее описания и скорости сходимости. Поэтому важность данной области возрастает со временем, оставаясь вне обсуждения в отечественной физикохимической литературе. Исключение составляет, вероятно, только метод эффективных потенциалов фрагментов (EFP) [7] вследствие популярности пакетов Firefly [8] и GAMESS [9], использующих EFP. Поэтому мы не будем рассматривать литературу, относящуюся к методу EFP. Отметим, что метод позволяет включать АММ до третьего порядка в декартовом представлении, которые можно рассчитать в рамках самого пакета.

Ниже мы рассмотрим методы распределенного мультипольного анализа (РМА) для молекулярных и твердотельных систем. Данные методы позволяют рассчитать атомные мультипольные моменты (АММ) нулевого порядка (заряды) и более высоких порядков.

1. МЕТОДЫ РАСПРЕДЕЛЕННОГО МУЛЬТИПОЛЬНОГО АНАЛИЗА

1.1. Методы распределенного мультипольного анализа для молекулярных систем

Модель точечного распределения электронной плотности (ЭПЛ) для описания электростатического потенциала (ЭПО) была впервые, по-видимому, реализована Халлом сначала для плавающих сферических гауссовых орбиталей s-типа на примере LiH [10, 11], а потом и для p-типа [12]. Идея представления молекул в виде характерных групп, переносимых между молекулами, для описания ЭПО была связана Томази и др. [13, 14] с переносимыми локализованными орбиталями (электронными парами), положения центров для которых строятся в иерархической последовательности до достижения требуемой точности описания ЭПО. В первых подходах центры для построения АММ были сдвинуты от атомных на расстояние до 3 а.е. [13–15], но позже Бентли показал достоинства схемы с атомными центрами [16]. Малликеновское распределение зарядов центров перекрывания (поровну между ближайшими атомами) было использовано в большей части указанных методов вместе с подходом Шипмана, который требовал сохранения дипольного молекулярного момента [17]. Затем Стоун использовал модель точечного распределения центров перекрывания [10, 11] и на базе теории сложения моментов [18] получил выражения для случая произвольной комбинации пар и типов (s‑, p‑, d-, f-) орбиталей, симметрии молекул с учетом более широкой выборки типов центров (и центров связей помимо атомных центров) [6].

Преимущества перехода от одно- к многоцентровой модели молекулы для представления электростатической энергии ее взаимодействия Uel связаны со значительным расширением области, в которой может быть получена аккуратная оценка Uel. Этот эффект иллюстрируется основным уравнением теории РМА для связи моментов [6]. Оно, например, может относиться к расчету центрального момента L-порядка изолированной молекулы $Q_{L}^{m}(A)$, заменяющего серию АММ на других (атомных) i-центрах молекулы:

(1)
$\begin{gathered} Q_{L}^{m}(A) = {{\sum\limits_{S = 0}^L {\sum\limits_{P = - S}^S {\left[ {\left( \begin{gathered} L + m \\ S + P \\ \end{gathered} \right)\left( \begin{gathered} L - m \\ S - P \\ \end{gathered} \right)} \right]} } }^{{1/2}}} \times \\ \times \;Q_{S}^{P}(i)R_{{L - S}}^{{m - P}}(A,i), \\ \end{gathered} $
где $R_{{L - S}}^{{m - P}}(A,i)$ = ${{\left( {\frac{{4\pi }}{{2(m - P) + 1}}} \right)}^{{1/2}}}{{r}^{{m - P}}}Y_{{L - S}}^{{m - P}}(\vartheta ,\varphi )$ – шаровые сферические гармоники, зависящие от вектора ${{\vec {d}}_{{iA}}}$ ($\left| {{{{\vec {d}}}_{{iA}}}} \right| = r$) между произвольным i-центром и атомом А. Это уравнение отвечает замене АММ $Q_{L}^{m}(A)$ порядка L серией АММ $Q_{S}^{p}(i)$ разных порядков S = 0, …, L на произвольном i-центре, таких, что их потенциал эквивалентен потенциалу, создаваемому одним $Q_{L}^{m}(A)$ [6]. Чем меньше расстояние r для переноса АММ в уравнении (1) относительно расстояния до точки, в которой оценивается потенциал, тем меньше характерное расстояние распределения заряда для данного центра и больше зона сходимости мультиполь-мультипольных взаимодействий. Хорошей иллюстрацией расширения зоны применимости мультипольных разложений является быстрый мультипольный (FMM) метод (см. ниже, часть 1.2.4) – постепенное дробление пар мультиполей, для которых достигнутое разделяющее расстояние между центрами ячеек является слишком коротким относительно размеров ячейки, на более мелкие зоны. Для мультиполей меньших ячеек достигнутое расстояние оказывается достаточно велико, что позволяет описать аналитически часть взаимодействий и тем уменьшить долю прямых КМ (квантово-механических)-расчетов.

Переход к атомным и связевым центрам связан с методом разделения пространства между выбранными центрами/атомами для расчета АММ. Вытекающее из него отнесение зарядовой плотности является одним из ключевых вопросов для всех методов, которые можно разделить на пространственно-ориентированные и базисно-ориентированные. К первой группе относятся варианты гиршфельдовской схемы распределения ЭПЛ (stockholder) [19–26], методы полиэдров Вороного [27], разделения по ван-дер-Ваальсовым (вдВ) атомным сферам [28], по поверхности промолекулы, отвечающей нейтральным атомам [29–31] и концепция поверхности “нулевого тока” (AIM, атом в молекуле) по Бейдеру [29–34]. Во вторую группу входят левдинская [35], малликеновская [36–39], стоуновская [6, 40], вейнхольдовская (NPA) [41, 42] схемы. С одной стороны, авторы, использующие интегрирование деформационной плотности, ставят под сомнение возможность базисно-ориентированных методов [27], особенно при переходе к плоским волнам [43]. Они аргументируют преимущества пространственно-ориентированных расчетов согласием относительных рассчитываемых зарядов со шкалой электроотрицательности по Полингу. При этом порядок рассчитываемых АММ ограничен нулевым [20–22, 27], первым [25, 26] или до квадруполей [23, 24], хотя для описания ЭПО, который является наблюдаемой величиной, с контролируемой точностью требуются мультиполи более высоких порядков для гиршфельдовской [44] малликеновской [1], стоуновской [44], КАММ (метод кумулятивных атомных мультипольных моментов) [44, 45] и бейдеровской [34] схем. Стоун обращал внимание на то, что асимметрия отдельных атомных фрагментов представляет наибольшую проблему для указанных подходов с интегрированием деформационной плотности, так как требует моментов высоких порядков [40]. Медленную сходимость ЭПО относительно увеличения порядка АММ для гиршфельдовской схемы отмечал Спакман [44], указывая, что и с добавлением гексадекаполей полученная величина хуже описывает полную энергию на КМ-уровне, соответствующем волновым функциям, использованным для РМА.

1.1.1. Пространственно-ориентированные методы РМА. Из методов пространственно-ориентированной группы контролируемая точность в расчете электростатического поля (ЭП) и ЭПО относительно требуемого верхнего порядка используемых АММ протестирована для разделения молекулы по поверхности “нулевого тока” электронной плотности (метод “атомов в молекуле”, АВМ или AIM) по Бейдеру, т.е., выбор отвечает поверхности, которую не пересекают линии тока. Поэтому результаты данного AIM-разделения пространства мы рассмотрим подробнее. Первая работа с использованием критериев AIM для расчета АММ была выполнена в 1985 г. [46], а впоследствии развивалась Ладигом [32], Попелье и колл. [29–31, 34] и многими другими. Работы первого направления были выполнены в декартовых координатах для моментов до третьего порядка, включая как электростатические моменты [32], так и тензоры поляризуемости [33]. Позже эти подходы были использованы при создании схемы ТАЕ (Transferable Atom Equivalents) для описания ЭПО вокруг относительно больших молекул, но схема показала только полуколичественные результаты [47]. В последней работе был предложен способ экономного интегрирования в рамках AIM-метода для произвольной структуры. Соответствующая МНК (метод наименьших квадратов) – ошибка составляет 25–36% в величине ЭПО, оцененного на поверхности постоянной электронной плотности вокруг молекул аланина, ди- и триглицина [47].

Также в декартовых координатах и близкой (гиршфельдовской) схеме распределения электронной плотности (ЭПЛ) была предложена альтернативная схема расчета АММ до второго порядка через интегрирование функции от электрического поля на поверхности “нулевого тока” [19], следуя идее Попелье рассчитывать так заряд вместо традиционного интегрирования негладкой функции ЭПЛ по объему атома (теорема Остроградского–Гаусса) [48]. Было выполнено сравнение диполей и квадруполей по поверхностям обоих типов (гиршфельдовской и бейдеровской), которое показало более аккуратные величины диполей с бейдеровской поверхностью и согласованный в обоих подходах существенный рост диполей (для HCN +33% и +31% по гиршфельдовской и бейдеровской поверхностям, соответственно) и квадруполей (для бензола +11.0% и +14.2% по гиршфельдовской и бейдеровской соответственно) всех рассмотренных молекул при переходе из газовой фазы в кристаллическую [19].

Активное продолжение анализ АММ получил с разделением ЭПЛ по поверхности “нулевого тока” в сферических координатах [29–31, 34]. Наиболее принципиальным моментом является высокая точность 0.1 кДж/моль для сходимости ЭПО с относительно невысокими порядками АММ, построенных по поверхности “нулевого тока” [34]. Такая сходимость требует учета АММ до третьего порядка для неводородных атомов и всего до второго порядка для водородов (по схеме только на атомных центрах) [34]. Аналогичная достаточность АММ до третьего порядка была показана ранее для гиршфельдовской, КАММ и стоуновской схем для сложного случая (HF)2 [44]. Необходимость увеличения числа центров относительно числа атомов была проанализирована ранее в соответствии с критерием для линейных молекул: Ω > 3Θ2/4μ, где μ, Θ, Ω обозначают дипольный, квадрупольный и октупольный моменты молекулы (z-компоненты) [15]. Усовершенствование схемы с включением центров на связях [30] не было широко применено в более поздних работах. Объяснение этому можно получить из работы [30]. Сравнивая РМА при малликеновском типе распределения ЭПЛ и РМА/AIM подходы для 10 вдВ-комплексов и 27 пар аминокислот, Жубер и Попелье обнаружили [30], что можно получить похожие по точности результаты, не используя для РМА/AIM-подхода точек на связях, хотя такие точки связей были учтены для традиционного РМА. Увеличение числа точек более одной на связях (АIМ1/4, АIМ1/6, АIМ1/6) или перенос их вне линий связей (модель АIМlap) привел к лишь немногим более качественному описанию энергии пар в рамках РМА/AIM [30]. Важной характеристикой для РМА подхода только по атомным центрам при разложении ∼ R–L до суммы взаимодействующих моментов шестого порядка (L = l1+ l2 + 1 = 6, где порядки АММ l1 и l2 относятся к разным молекулам) были средние ошибки в энергии около 0.24 ккал/моль для набора вдВ кластеров и 1.2–1.3 ккал/моль для пар аминокислот [30]. Интересно, что худшие результаты РМА/AIM модели по всем параметрам сходимости получены в случае учета одной дополнительной точки на связях (АIМ1/2) по сравнению с только атомными центрами (АIМ0) [30]. В традиционной РМА схеме ошибки уменьшаются, а сходимоcть растет при расширении числа центров [30]. Комментируя данную работу [30], Стоун отмечал, что в конденсированной фазе ЭПЛ благодаря соседним молекулам более сжата, а потому АММ более локализованы [40]. Тем не менее, это замечание не объясняет эффект ухудшения сходимости РМА/AIM-схемы [30] при включении центров разложения на связях. Одновременный расчет авторами стоуновских АММ с помощью пакета GDMA (распределенный мультипольный анализ) [49] и их анализ показал ожидаемое улучшение по точности и скорости при введении одной связевой точки для серии маленьких систем (вдВ-кластеров). Для больших систем (пар аминокислот) сходимость достигается и без новых точек, приводящих к меньшей точности. Разложения по стоуновским АММ для больших систем оказались и более аккуратными и быстрыми [30].

В других работах группы Попелье и сотр. [29] исследуется возможность построения “топологических” силовых полей (СП), соответствующих описанию кулоновской ЭЭВ (электростатическая энергия взаимодействия) в рамках распределения ЭПЛ по Бэйдеру, т.е., атомы в молекуле (АВМ или AIM). Показано, что в одной молекуле при увеличении порядка АММ до 4, сходятся с разумной точностью (до 26% в бутане для пары С1–С3) взаимодействия атомов, разделенных хотя бы одним атомом при исключении цепей содержащих кратные связи (акролеин [29]). Взаимодействия ближайших соседей в одной молекуле не сходятся при увеличении порядка АММ до 4 [29], а между атомами разных молекул – сходятся. Чтобы достичь количественной сходимости соседей 1–3- и 1–4-типов в одной молекуле, была выполнена попытка выразить те же взаимодействия через моменты, смещенные по уравнению (1) [6] так, чтобы увеличить расстояние между новыми центрами разложения ЭПЛ и ускорить сходимость [30]. Алгоритм предполагает сдвиг центров в противоположных направлениях относительно друг друга. Несмотря на необходимость моментов более высоких порядков, такой подход приводит к успеху.

Новым направлением развития моментного описания в рамках АIМ стало применение Попелье и сотр. схемы непрерывных AММ, полученных в целых бесселевых функциях (нулевого порядка) J0(r) [31], например, момент атома А с областью:

(2)
$Q(t) = \int\limits_{{{\Omega }_{A}}} {dr{\kern 1pt} {\text{'}}\rho (r{\kern 1pt} '){{e}^{{r't\cos \gamma }}}{{J}_{0}}(r{\kern 1pt} {\text{'}}t\sin \gamma ),} $
где γ – угол меду вектором r, задающем точку расчета ЭПО, и вектором r' по зарядовому распределению атома А, а ЭПО, создаваемый атомом А, при этом численно интегрируется с весами Гаусса–Лежандра wi на сетке из nquad точек в виде:

(3)
$V{{(r,{{\Omega }_{A}})}_{{elec}}} = \sum\limits_{i = 1}^{{{n}_{{quad}}}} {{{w}_{i}}{{e}^{{ - r{{t}_{i}}}}}Q({{t}_{i}})} .$

Параметр t предполагается изменяющимся в конечных пределах [0, T] вместо интервала [0, ∞], формально получающегося при разложении |r – r'|–1, что позволяет оценивать потенциал при умеренном порядке числа моментов и меры Т. Для каждого направления расчета ЭПО моменты должны рассчитываться заново, их порядок определяется требуемой точностью описания ЭПО. Tо есть, заранее следует определить число выбранных направлений, в которых будет рассчитан ЭПО, и для каждого направления должны быть получены наборы АММ. Для сравнения, при nquad = 20 и Т = = 5 в представлении ЭПО (3), относительная точность описания ЭПО, создаваемого атомом N в молекуле N2 или атомом О в молекуле Н2О, достаточно высока и составляет 10–2 кДж/моль на границе области атома в АIМ модели (около 3.4 Å от центра N2 молекулы) [31].

Гибридным вариантом между базисно-ориентированными и пространственно-ориентированными методами РМА стали подходы на основе ISA [21, 22], представленные в двух вариантах – Minimal Basis Iterative Stockholder (MBIS) [50] и basis space ISA (BS-ISA) [51]. В обоих случаях полная электронная плотность на атоме ρa(r) определяется по всем волновым функциям, а вес атома wa(r) в полной ЭПЛ определяется по сферическим функциям (см. ниже, часть 1.1.4).

1.1.2. Метод РМА по Стоуну. Наиболее распространенной схемой РМА из базисно-ориентированных методов является схема Стоуна с отнесением зарядовой плотности к ближайшему соседему атомному центру (или связевому, если определен) [6], чем отличается от малликеновской схемы (разделение электроннной плотности пополам между двумя соседними центрами) [36–39]. Одним из основных проблемных мест данной схемы, как и других базисно-ориентированных, является зависимость АММ от выбора базиса. Качественные изменения (на уровне геометрии молекул) в электростатических взаимодействиях при изменении базиса были отмечены [52] много позже появления первых результатов, указывающих на сильную зависимость величин АММ от базисов [45, 53–55]. В системах с водородными связями корректное описание квадруполей кислородов приводит к точному описанию геометрии [52]. Для схемы Стоуна изменение ЭЭВ в H2CO…HF кластере было рассмотрено в ряду пяти базисов STO-3G, 3-21G, 6-31G, 3-21G* и 6-31G* [52]. Требуемый для согласия с экспериментом двухъямный потенциал ЭЭВ относительно угла поворота HF по отношению к плоскости H2CO был получен, только начиная с 6-31G базиса. Аналогично роль квадруполей атомов углерода в димере бензола определяет Т-образную геометрию, наблюдаемую в кристаллах [52].

Это важное достоинство данной схемы РМА – способности описывать качественные особенности структуры молекулярных кристаллов – распространяется на широкий диапазон объектов. Так электростатическое взаимодействие на малых расстояниях при упаковке кристаллов в зависимости от максимальных величин диполей или квадруполей, определяют пространственную группу симметрии кристаллов молекул Х2 [56]. Молекулы, принципиально отличающиеся по величинам АММ, кристаллизуются в разных группах пространственной симметрии. Для молекул галогенов, обладающих значительными атомными квадруполями вследствие неподеленных электронных пар, характерна группа Сmca, а для N2, имеющей центральный положительный квадруполь от электронов π-связей и существеные атомные диполи, кубическая группа Ра3 наиболее приемлема [56]. Для широкой группы полярных и полярных ароматических молекул использование перемасштабированных (*0.9) АММ изолированной молекулы на ХФ/6-31G** уровне, позволило получить оценки параметров соответствующих кристаллов, с учетом дисперсионной ЭВ (ДЭВ) в виде потенциала Бэкингэма [57]. Расширенное обсуждение моделирования органических кристаллов и других объектов с помощью уже разработанных к 2005 г. программ, DL_POLY и DL_multi оперирующих с АММ, можно найти в работе [58], анализирующей их применение в молекулярнодинамических расчетах.

Систематическое исследование АММ базиснозависимых методов мы можем найти в нескольких работах, выполненных в рамках схем Стоуна [49, 54] и КАММ [45]. В первом случае [54] показаны качественные изменения и АММ, и связевых ММ в молекуле N2 между базисом А (5s3p1d с диффузными 1s1p2d-функциями на центре N–N-связи), с одной стороны, и базисами B (5s4p3d и две поляризационные функции на центре N–N-связи) и C (как и B, но другой выбор поляризационных функций), с другой стороны. Между базисами B и C различия носят количественный характер. Для базиса B были найдены также только количественные изменения АММ на трех разных уровнях теории: ССП (метод самосогласованного поля), CAS (метод полного активного пространства) и МР2. Эти же три метода в применении к молекуле воды с базисом [6s5p2d/4s3p] показали качественные (атом Н) и количественные (связи и атом О) изменения моментов до квадруполя во всех положениях, хотя сами диполи и квадруполи на связях и атоме Н невелики по величине [54].

Система РМА по Стоуну, менее зависимая от смены базиса, была построена Уитли для линейных молекул [59]. Возможность выбора такой системы центров была ранее высказана Стоуном [6]. Выражения для мультиполей более высоких порядков $Q_{{l'}}^{{{\text{new}}}}({{z}_{p}})$, ограниченных по порядку l' сверху как lmax > l' и распределенных по р-точкам, p = 1, …, n, через известный момент Ql(zA) более низкого порядка (l' ≥ l), например, молекулярный, отнесенный к центру zA молекулы (к центру заряда для иона), были получены [59] в рамках РМА схемы Стоуна [6]:

(4)
$\begin{gathered} Q_{{l'}}^{{new}}({{z}_{p}}) = \left( \begin{gathered} l{\kern 1pt} ' \\ l \\ \end{gathered} \right)\frac{{({{l}_{{\max }}} - l){\kern 1pt} !({{l}_{{\max }}} - l{\kern 1pt} {\text{'}} + n - 1){\kern 1pt} !}}{{({{l}_{{\max }}} - l{\kern 1pt} {\text{'}}{\kern 1pt} ){\kern 1pt} !({{l}_{{\max }}} - l + n - 1){\kern 1pt} !}} \times \\ \times \;{{({{z}_{{\text{A}}}} - {{z}_{p}})}^{{l' - l}}}{{Q}_{l}}({{z}_{{\text{A}}}})\frac{{{{\Pi }_{{p' \ne p}}}({{z}_{{\text{A}}}} - {{z}_{{p'}}})}}{{{{\Pi }_{{p' \ne p}}}({{z}_{p}} - {{z}_{{p'}}})}}. \\ \end{gathered} $
Для этого Уитли была доказана лемма о поведении соответствующей функции в выражении (4). Из $Q_{{l'}}^{{{\text{new}}}}$ рассчитывались центральные моменты системы $Q_{L}^{{{\text{tot}}}}$:
(5)
$Q_{L}^{{{\text{tot}}}}({{z}_{{\text{A}}}}) = \sum\limits_p {\sum\limits_{l' < L} {\left( \begin{gathered} L \ \\ l{\kern 1pt} ' \ \\ \end{gathered} \right)} {{{({{z}_{p}} - {{z}_{{\text{A}}}})}}^{{L = l'}}}Q_{{l'}}^{{{\text{new}}}}({{z}_{p}}).} $
Учет отброшенных мультиполей более высоких порядков l > lmax допускался за счет модификации только максимально разрешенного $Q_{{l\,\max }}^{{{\text{new}}}}$ в (4) [59]. Данный метод позволил получать менее зависимые от базиса величины центральных моментов из распределенных $Q_{{l'}}^{{{\text{new}}}}$, хотя сами распределенные АММ могли существенно меняться, например, от 8 до 172% для АММ первых (L = 0–2) низших порядков на Н-атоме в системе (HF)2 и от 8 до 87% на F-атоме.

Выражение (5) наиболее близко из представленных в литературе к предложенному одним из авторов статьи (АВЛ) выражению для поиска корреляции зарядов и АММ:

(6)
$Q_{L}^{m}(A) = {{a}_{L}}_{m}R_{L}^{m}(A) + {{b}_{{Lm}}},$
где aL и bL рассматривались как подгоночные параметры, а функция $R_{L}^{m}(A)$ выражается через заряды (в отличие от (5), где используются все моменты низших порядков l' < L) всех соседей атома А:
(7)
$R_{L}^{m}(A) = \sum\limits_{i = 1}^N {Q_{0}^{0}(i)X_{L}^{m}(A,i)d_{{i{\text{A}}}}^{{ - K}}} ,$
где $X_{L}^{m}(A,i)$ соответствует ненормированным полиномам Лежандра (выражения приведены в [60]), применяемым в рамках метода [1], K = GL + + 1, ${{\vec {d}}_{{iA}}}$ – вектор между произвольным i-центром и атомом А, ${{d}_{{i{\text{A}}}}}$ = $({{({{X}_{i}} - {{X}_{{\text{A}}}})}^{{\text{2}}}} + {{({{Y}_{i}} - {{Y}_{{\text{A}}}})}^{{\text{2}}}} + ({{Z}_{i}} - $ $ - \;{{Z}_{{\text{A}}}}{{)}^{{\text{2}}}}{{)}^{{{\text{1}}/{\text{2}}}}}$. Используя выражения (6), (7) которого удалось получить высокую корреляцию “$R_{L}^{m} - Q_{L}^{m}$” для представления АММ $Q_{L}^{m}$ полностью кремниевых цеолитов [61], включая аморфные [62], алюмофосфатов [63–65] и катионных форм [66]. Отношение позволяет экстраполировать АММ на произвольную структуру того же состава по известной геометрии. Мы описали данные результаты более подробно ниже в части 2.

Частичное решение для произвольных молекул по преодолению зависимости РМА по Стоуну от выбранного базиса, предложенное в виде интегрирования на сетке [49], аналогичное сетке для DFT (метод теории функционала плотности)-интегрирования, но в более гибкой форме отнесения зарядов – в виде перекрывающихся (fuzzing) полиэдров Вороного [67] с частичным учетом симметрии [68]. Модифицированная схема Стоуна реализована в виде пакета GDMA [49]. В то же время это усовершенствование требует значительного времени счета, в первую очередь, для диффузных функций.

Многие из проблем переносимости АММ между разными молекулами пептидов в методе Стоуна обсуждалась в работах Прайс и колл. [53, 69, 70]. Авторы предположили, что первая проблема – изменение электронного состояния отдельных групп (как O=C–N) и зарядов вследствие изменения конформации (например, вследствие перераспределения зарядов по механизму гиперконьюгации в плоской модели дипептида) – может потребовать в дальнейшем эмпирических поправок, то есть, выхода за рамки обычного РМА. Решение второй важной для РМА проблемы – учета поляризуемости групп с изменением конформации – было предложено и протестировано с помощью распределенных поляризуемостей [53]. Авторы [53, 70] предполагали, что не водородная связь (ВС) является основной проблемой для переносимости, а зависимость АММ от 1–4 торсионных углов. Применимость АММ сравнивалась между тремя конформерами N-ацетилаланина-N'-метиламида (NAAla) и четырьмя конформерами N-ацетилдиглицида-N'-метиламида (NAdGlu) [53]. Авторы [53] критиковали заключение о переносимости зарядов из монопептидов на полипептиды, данное в работе Пульмана и колл. [71] из-за ограниченных отличий между тестовыми фрагментами (Lys, Asp, Tyr, Gln) и тренировочной молекулой. Они пришли к выводу о приемлемой точности предположения о переносимости АММ, полученных из дипептидов [53, 69] в согласии с работой [72]. Работая с дипептидом из одинаковых исходных аминокислот (диглицид), авторы не рассматривали усложнение в обычных смешанных полипептидах, связанное с разделением заряда между звеньями разных кислотных остатков и соответствующий этому возможный сдвиг всех моментов, считая такое влияние несущественным [53] или нейтрализуя разницу зарядов при переносе из кристаллической фазы без обсуждения деталей (меняя заряд только атома Сα между пептидными остатками), как в их более ранней работе [69]. Там были оценены ненулевые заряды фрагментов в дипептидах, но их влияние было оценено как крайне слабое. Как существенная была выделена лишь “основная” ошибка, быстро падающая с расстоянием от атомов, при расчете ЭПО без учета изменения торсионных углов (относительно геометрии дипептидов в кристаллической фазе) возле протона амино-группы с ошибкой 7–27 кДж/моль и чуть меньшей у карбонильного кислорода (9–17 кДж/моль, обе точки на расстоянии 2.06 Å от атомов). Эта ошибка от эффекта поляризации атомов основной цепи боковой цепью возле атома азота снижается от 0–10 кДж/моль (3.44 Å) до 0–3 кДж/моль (8.44 Å). В этой работе поляризация атомов основной цепи атомами боковой цепи рассматривалась с базисом 3-21G без поляризационных функций, поэтому эффект мог быть занижен [69]. Но более ионная структура с 3-21G (средние заряды $Q_{0}^{0}({\text{C}})$ = 1.050 e, $Q_{0}^{0}({\text{O}})$ = –0.917 e по группе дипептидов с боковыми ветвями, образующими спираль, табл. 1 из [69]), по-видимому, меньше влияет на ЭПО на больших расстояниях от молекулы. Для сравнения, применение 3-21G и 6-31G** привело к сравнимым ошибкам на поверхности ЭПО около NAAla или NAdGlu, существенно удаленной от молекулы (1.4 Å над ван-дер-Ваальсовой поверхностью) [53]. Обычно приводящий к высоким малликеновским зарядам базис 3-21G может быть причиной того, что заряд (±1) был полностью сосредоточен на боковой цепи и не влиял на заряды атомов основной цепи. Учет заряда монопептидного фрагмента (остатка) –NHRCHRCO – позже был изучен Сокальским и колл., которые показали важность сохранения корректного зарядового баланса на отдельных остатках, если АММ и заряды переносятся из произвольной библиотеки пептидов [73] (см. ниже). Заряды, рассчитанные на ХФ и B3LYP уровнях с базисом 6-31G**, которые распределяли (по остаткам) Сокальский и колл. (–0.0463, –0.0358 е на Ala и Gly соответственно), были того же порядка, что и в работе [70]: +0.046, +0.075, и +0.042 e (базис 3-21G) на N-метилированных остатках Gly, Val и Leu, где последствия нейтрализации зарядов остатков не оценивались. Возможно, что различие оценок влияния остаточных зарядов фрагментов является следствием разных схем оценки АММ по Стоуну и в КАММ.

Таблица 1.  

Обозначения и типы использованных базисов для расчета ЭПО и ЭП в цеолите MgPHI (цеолит филлипсит) (Т = Si, Al)

Номер базиса Базис
BS1 STO-3G
BS2 86-1G(Mg)/6-21G(T)°/6-311G(O)
BS3 86-1G(Mg)/85-11G(Al)/6-21G(Si)/8-411G(O)
BS4 86-1G(Mg)/6-21G(T)/6-311G*(O)
BS5 85-11G*(Mg, Al)/8-31G*(Si)/6-31G*(O)
BS6 85-11G*(Mg, Al)/6-21G*(Si)/8-411G*(O)
BS7 85-11G(Mg, Al)/8-31G(Si)/8-411G*(O)
BS8 85-11G*(Mg)/8-31G*(T)/6-31G*(O)
BS9 86-1G*(Mg)/6-21G*(T)/6-311G*(O)
BS10 85-11G*(Mg, Al)/8-31G*(Si)/8-411G*(O)
BS11 86-1G*(Mg)/8-31G*(T)/6-311G*(O)
BS12 85-11G*(Mg)/8-31G*(T)/6-311G*(O)
BS13 86-1G*(Mg)/85-11G*(Al)/8-31G*(Si)/8-411G*(O)
BS14 85-11G*(Mg)/8-31G*(T)/8-411G*(O)

1.1.3. Метод КАММ. В рамках последнего метода (КАММ) Сокальского [45], оперирующего с декартовыми координатами и только атомными центрами (без центров связей), знание элементов матрицы плотности Prs, описывающей связи с соседними атомами, необходимо для последовательного построения АММ $m_{c}^{{klm}}$ из зарядов Zс и АММ $m_{c}^{{k'l'm'}}$ низших порядков того же самого с-атома:

(8)
$\begin{gathered} m_{c}^{{klm}} = {{\delta }_{{ac}}}{{Z}_{a}}u_{c}^{k}{v}_{c}^{l}w_{c}^{m}\sum\limits_{r \in c} {\sum\limits_s {{{P}_{{rs}}}\left\langle {r\left| {u_{c}^{k}{v}_{c}^{l}w_{c}^{m}} \right|s} \right\rangle } } - \\ - \;\sum\limits_{k' > 0} {\sum\limits_{l' > 0} {\sum\limits_{m' > 0,(k',l',m' \ne k,l,m)} {\left( \begin{gathered} k \ \\ k{\kern 1pt} ' \ \\ \end{gathered} \right)\left( \begin{gathered} l \ \\ l{\kern 1pt} ' \ \\ \end{gathered} \right)\left( \begin{gathered} m \ \\ m{\kern 1pt} ' \ \\ \end{gathered} \right)} } } \times \\ \times \;u_{c}^{{k - k'}}{v}_{c}^{{l - l'}}w_{c}^{{m - m'}}m_{c}^{{k'l'm'}}, \\ \end{gathered} $
где декартовы координаты $(u,{v},w)$ = $(x,y,z)$, $\left( \begin{gathered} k \ \\ k{\kern 1pt} ' \ \\ \end{gathered} \right)$ – биномиальные коэффициенты. Из выражения (8) вытекает “кумулятивность” подхода – все АММ данного атома высокого порядка строятся из его же АММ более низких порядков. Как и в случае РМА Стоуна, при неизвестных величинах Prs в (8) или $P_{{{\text{12}}}}^{{\text{g}}}$ в выражении (12) ниже, которые можно получить только из прямого расчета, нельзя экстраполировать величину КАММ для аналогичной системы. Данный подход не ограничен одним способом разнесения электронной плотности по атомам, и необходимость использования помимо зарядов и КАММ высших порядков для количественной оценки ЭПО была показана как с левдинскими [35] и малликеновскими [36–39] зарядами, так и с зарядами, оцененными из ЭПО (PDC (метод расчета зарядов, подогнанных по электростатическому потенциалу) или ESP (метод расчета зарядов, подогнанных по электростатическому потенциалу) типов зарядов) [73, 74].

Как базиснозависимый метод, КАММ при переходе между базисами 6-31G* и 6-311G(2d/2p) для молекулы HCN с малликеновской схемой привел к резкому изменению атомных зарядов и диполей (со сменой заряда компоненты $m_{С}^{{00{\text{1}}}}$ диполя на атоме С), хотя молекулярный диполь HCN изменился несущественно [74]. В случае СО, крайне чувствительной к учету корреляции, с малликеновской схемой переход между базисами 6-31G* и 6-311G(3d) существенно изменяет все АММ, но поскольку расчет был выполнен на уровне ССП ХФ, то диполь СО (–0.130 а.е. → → ‒0.088 а.е.) так и не приближается к экспериментальной величине (0.057 а.е. [75]) [74]. Для HCN было проверено поведение при изменении ССП ХФ уровня расчета на КВ (с однократными и двукратными возбуждениями или CISD), что привело к умеренным изменениям в величинах КАММ и улучшению согласия с экспериментальной величиной диполя HCN. При смене ХФ → → CISD с базисом 6-311G(2d/2p) изменение составило 6.8%, что близко к аналогичному сдвигу при другом выборе зарядов для схемы КАММ (подгонка зарядов к ЭПО) и базисе 6-31G* (6.4%) [74].

Наиболее полные наборы базисов (до шести типов для N2 и СО) были изучены в работе для оценки вкладов электронной корреляции (на уровне CISD, CAS, MP3, MP4SDTQ) в молекулярные моменты (до октуполя), рассчитанные через КАММ. Сами КАММ в явном виде приводились только для LiH и вместе с графическими результатами обрезания КАММ по разным порядкам для расчета ЭПО на оси С2Н2 [45]. При учете электронной корреляции на уровне конфигурационного взаимодействия с одно- и двухкратными возбуждениями (CISD) при базисе 6-311G(3d,2p) изменения КАММ в LiH не превышали 10% относительно некорреллированных величин и были того же порядка, что и изменения молекулярных моментов. Например, для LiH по абсолютной величине КАММ изменения составили –8.5% для заряда (m000), 1.4% (Li), –6.7% (H), –2.3% (LiH) для компонент m001 диполя, 2.5% (Li), –5.2% (H), 9.1% (LiH) для компоненты m200 квадруполя.

До начала обсуждения в литературе зависимости АММ (и КАММ) от конформационных углов [52, 53, 69, 70] автором [76] была составлена библиотека КАММ (до квадруполя) для ряда 28 молекул аминокислот. В более позднем приложении метода КАММ к циклическим производным аминокислот аланина и глицина типа (Ala)6 – x(Gly)x, x = 0, 2, 4, была тестирована переносимость АММ (до октуполя), рассчитанных на тестовых молекулах двух типов – либо аналогичных циклических производных, либо более широкой библиотеки производных аминокислот (Glu, His, Lys и Arg), включая заряженные фрагменты [73]. Между похожими циклическими производными (Ala)6 – x(Gly)x (различие в зарядах отдельных кислотных фрагментов в тестовой и исследуемой системах не превышало 0.0068 е) ошибка в ЭПО при перенесении АММ уменьшалась до 6% (с учетом членов до показателя R–4 в энергии взаимодействия) для разных типов циклических молекул (x = 0, 2, 4) с корректировкой до квадруполя включительно. Скорее неудачей окончилась попытка улучшить ОМНК (метод относительных наименьших квадратов)-разницу (уже выше показателя R–1 в энергии взаимодействия) с ошибкой одного порядка около 16–42% (различие в зарядах 0.036–0.046 е отдельных кислотных остатков Ala, Gly) при рассмотрении тестовых молекул из библитотеки аминокислот. Четыре предложенных схемы коррекции зарядов привели к уменьшению ошибки только на уровне взаимодействия заряд–заряд, но более высокие члены энергии взаимодействия не удалось описать по КАММ схеме при коррекции только зарядов, а из них – АММ по выражению (1.3). Такое различие зарядов даже атомов (фрагменты указаны в [73]) одного элемента (во втором знаке), как во втором случае работы [73], между разными оптимизированными конформерами дипептида аланина является вполне типичным, как следует из результатов рассмотрения в рамках метода Стоуна [70]. Это указывает, что проблема согласования зарядов пептидных остатков для переносимости АММ между пептидами будет возникать при использовании любых библиотек данных.

1.1.4. Методы РМА с ограниченным порядком АММ. Переопределенность полной системы АММ хорошо известна [77] и представляет одну из проблем схем РМА для расчета ЭП и ЭПО, для решения которой предлагаются возможные способы обрезания АММ высших порядков и сохранения наиболее важных АММ низших порядков. Один из подходов к выбору минимального числа АММ (вместе с зарядами), достаточного для описания ЭПО возле изолированной молекулы, является чисто комбинаторным [77, 78]. Делаются оценки доли включаемых в описание АММ от их полного набора [78], достаточной для относительно малой ошибки описания ЭП, напр., 0.7 ккал/е, вокруг молекул аминокислот из 19–25 атомов [77]. Текущей задачей можно считать поиск алгоритмов, которые позволили бы перейти к менее чем 2N моментам (вместе с зарядами) при N атомах в молекуле [78].

Число используемых центров может существенно превышать N в методах, ориентированных на подгонку только АММ нулевого порядка из условия описания ЭПО – зарядов (ESP или electrostatic potential charges) – на достаточно плотной сетке рядом с изучаемой системой. Если число требуемых зарядов для воспроизведения ЭПО ограничено только числом N атомных центров, то формально такая система зарядов может использоваться и для экстраполяции для новой системы. При больших N силовым полем CHELPG атомные заряды подгоняются на сетке точек, расположенной между поверхностями с общими атомным центром, ограниченными атомным ван-дер-Ваальсовым радиусом и радиусом 2.8 Å [79] или, по крайней мере, радиус сеток должен превышать ван-дер-Ваальсов (вдВ) радиус в 1.2 раза. Для тестирования могут сравниваться величины ЭПО и на сферах с превышением вдВ) радиуса в 1.5 и 2 раза [51]. Увеличение числа точек на порядок в последней работе [79] по сравнению с изначально предложенной сеткой в CHELP модели [80, 81] позволило обеспечить гладкость рассчитываемого потенциала относительно конформационных поворотов молекулы. Из-за большого числа требующихся точек более экономным оказался способ подгонки атомных зарядов к величинам, позволяющим описать молекулярные мультипольные моменты до октуполя [82] или квадруполя (только с атомными центрами для зарядов [83] или сеткой зарядов возле каждой атомной позиции [84]), а не ЭП. Наиболее распространенной техникой для подгонки зарядов и АММ стал метод множителей Лагранжа, начиная с применения для подгонки зарядов в моделях ESP [85], CHELP [80, 81], CHELPG [79], RESP [86] в более поздней модели REPEAT [87], а также для воспроизведения ЭПО или молекулярных мультипольных моментов [82, 83]. Недостаточность представления только центральных мультипольных моментов для воспроизведения поля и потенциала отмечалась Шпакманом [44]. Тем не менее, такой MDC (заряды, полученные из мультипольных моментов)-подход оказался удовлетворительным в работе [83]. Гиллом и др. [83] была использована идея аппроксимации молекулярных мультипольных моментов (МММ) зарядами, чтобы получить минимальный набор МММ, подгоняя их разложениями по атомным зарядам и гармоникам соответствующего порядка. Данный набор зарядов, подогнанных под описание мультипольных моментов (общепринятое обозначение MDC для этого типа), должен описывать как можно больше компонент МММ. Число зарядов может быть равно или меньше числа описываемых МММ. В рамках данной модели авторы показали более аккуратное описание ЭПО по сравнению с рядом известных способов представления зарядов (CHELPG, NPA, AIM, RESP, MK, Малликеновские заряды и др.) над и под плоскостью симметрии метанола на уровне B3LYP/6-31G(d) и в плоскости нитробензола (CHELPG и малликеновские заряды) на уровне ХФ/6-31G(d) вместо иллюстрации величин только МММ [82]. Преимущество данного метода по сравнению с натуральными зарядами (NPA) [41, 42] и малликеновскими зарядами было более очевидно для нитробензола. Но в этом случае физический смысл зарядов утерян, а, именно, заряды, оцененные из мультиполей (MDC) атомов N, О, С(–N), C (в пара-положении относительно ‒NO2-группы), оказались равны –0.800, 0.018, 3.247, 2.799 е соответственно. Анализ зарядов (малликеновские, гиршфельдовские, NPA, GAPT, CHELPG, AIM и др.) с акцентом на атомы Н и С на хартри-фоковском уровне выполнялся Вибергом и Рабленом [88], а использование зарядов в рамках разных схем и АММ низших порядов (до квадруполя) для построения силовых полей подробно описывалось Сиплаком [89].

Вместо метода множителей Лагранжа Ципер и Бурке [90] предложили свой подход минимизации числа атомных мультипольных моментов изолированной молекулы. Им стало решение системы линейных уравнений методом minimal atomic multipole expansion [90] с учетом мультиполей (до квадруполей) тех неподеленных атомных пар, которые вытекают из льюисовской модели молекулы. Метод основан на выделении вокруг молекулы изоповерхности, соответствующей достаточно большой доле ЭПЛ, зарядом вне которой можно было бы пренебречь. Тогда потенциал, создаваемый на данной поверхности, считается созданным только мультиполями, находящимися внутри данной поверхности. Из соответствующей линейной системы уравнений, отвечающих минимуму разности точного и приближенного ЭП на изоэлектронной поверхности, определяются АММ и заряды, находящиеся внутри данной поверхности.

Дополнительные центры связей (помимо атомных центров) были введены только для небольших молекул, таких, как H2О, Н2, HF, F2, CO [82]. В последней работе на наборе из 31 молекулы и нескольких пептидов, рассчитанных на DFT/BP/TZ2P-уровне, в функционал, зависящий от разности подгоняемых и рассчитанных АММ (отдельный множитель Лагранжа для каждого порядка до октуполей), был введен экспоненциальный множитель (wA,s), демпфирующий вклады удаленных соседних атомов (по сравнению с ближайшим атомом) и “локализующий” АММ. Этим масштабированием и формой функционала (от разности АММ) работа [82] и отличается от похожего (но много более раннего) подгонки АММ (до квадруполей включительно) для набора из 24 молекул на уровне ХФ/6-31G**. В этом случае АММ были введены в оптимизируемый функционал от разности ЭПО [91]. Автор представил доскональное сравнение семи вариантов подгонки от только зарядов до всех АММ вместе (заряды, диполи, квадруполи). Акцент сделан на связевые диполи, потому что уже на уровне “заряды плюс диполи” достигается ОМНК менее 3%. Обе системы диполей, направленных строго по связям или без ограничений, привели к хорошему описанию ЭПО молекул. Более подробно влияние дополнительных центров для центров связей и их положения обсуждалось для РМА и РМА/AIM подходов Жубером и Попелье (см. ниже, [30]). Возможности повышения качества представления ЭП с использованием малликеновских зарядов при условии расчета заряда по отдельным связям атома была показана в рамках расщепленного метода выравнивания электроотрицательностей [92]. Применение данного расщепленного (по связям) метода было более полно охарактеризовано для цеолитов с применением Гиршфельдовских зарядов [93].

Одна из проблем применения Гиршфельдовской концепции зарядов к неионным системам в том, что для ряда атомов нельзя описать стандартное состояние в рамках только одного электронного состояния [94]. Другую задачу – произвольный выбор начального стандарта (проатомов в промолекуле) – было предложено преодолеть итерационным путем в работе Бутлинка и др. [20], которые использовали целевую функцию (GF) Парра и др. [95], взятую из теории информации:

(9)
$GF = \int {{{\rho }_{{\text{A}}}}(r)\ln \left( {\frac{{{{\rho }_{{\text{A}}}}(r)}}{{\rho _{{\text{A}}}^{0}(r)}}} \right)} {\kern 1pt} dr + \int {{{\rho }_{{\text{B}}}}(r)\ln \left( {\frac{{{{\rho }_{{\text{B}}}}(r)}}{{\rho _{{\text{B}}}^{0}(r)}}} \right){\kern 1pt} dr} $
для образования молекулы АВ, и ввели функцию формы (или нормированной плотности σA= = ρA/NA) для строгой связи с выражением по теории информации. Близкий к (9) функционал GF был получен Уитли и Лиллестолен [21, 22] совместным решением уравнений для сферически-усредненной плотности атома 〈ρa(r)〉 и его веса wa(r) в общей электронной плотности ρ(r):

(10)
${{w}_{a}}({\mathbf{r}}) = \langle {{r}_{a}}({\mathbf{r}})\rangle ,$
(11)
${{\rho }_{a}}({\mathbf{r}}) = \rho ({\mathbf{r}}) \times ({{w}_{a}}({\mathbf{r}}){\text{/}} > {{\Sigma }_{i}}{{w}_{i}}({\mathbf{r}})).$

Такое решение (Iterative Stockholder Approximation или GISA) оказалось независимым от выбора стандарта для проатомов и получило наибольшее распространение. Оно применено как тест в методе эффективных потенциалов фрагментов (EFP) [96] в пакете CRYSTAL14 [4], послужило основой для методов Basis Set-ISA (BS-ISA) [51], Gaussian Iterative Stockholder Approximation (GISA) [96], Minimal Basis Iterative Stockholder (MBIS) [97], метода ISA-Pol для локальных поляризуемостей [98] и др. ISA вариант Уитли и Лиллестолен [21, 22] (вместе с традиционным Гиршфельдеровским) был использован Манцем и др. [23, 24], учитывая АММ до квадруполей, для расчета по наиболее широкому набору объектов из приводимых в литературе, содержащего твердотельные и пористые объекты.

Применение итерационной схемы [20] к молекулам и радикалам с ненулевым спином было показано в более поздней FOHI модели [26] в рамках Гиршфельдовского распределения ЭПЛ и GF (9). Работа [26] расширила число моментов до включения диполей и была протестирована сравнением на уровне теории UB3LYP, UMP2, UCCSD (базис aug-cc-pVTZ для UMP2 и UCCSD и 6-31G, 6-311++G** для UB3LYP) и с двумя эскпериментальными распределениями спиновой плотности. Был сделан вывод о корректной, но преувеличенной спиновой заселенности с FOHI-моделью. В ее поздней FOHI-D модификации [25] расчет зарядов и диполей выполнялся независимо. При этом последняя часть процедуры (расчет диполей) оказалась базисно-зависимой, что снижает преимущества Гиршфельдовского распределения. Авторы не сравнили новый вариант FOHI-D для случая с ненулевым спином [25]. Переносимость FOHI и FOHI-D атомных зарядов между системами аналогичного состава не была исследована.

1.1.5. Методы РМА из ЭПЛ на базе полуэмпирических методов. Разработка моделей АММ по результатам расчета полуэмпирической ЭПЛ велась с момента появления полуэмпирических методов.

Первые схемы РМА были применены к ЭПЛ, рассчитанной полуэмпирическими методами, такими как расширенный метод Хюккеля (IEHT) [99] или CNDO [99, 100]. В последней работе рассматривалась схема РМА для произвольной ЭПЛ с акцентом на слейтеровские АО в декартовом представлении, как и последующая схема, использованная авторами в пакете CRYSTAL [1] (см. часть 1.2.1), а применена была к ЭПЛ, рассчитанной на CNDO/2-уровне [100]. Анализировались четыре варианта распределения ЭПЛ, в том числе, с выделением АММ валентных электронов и учетом внутренних оболочек в виде зарядов, отнесенных к положению ядер. Такая схема оказалась самой неудачной из четырех, что отражает особенность ЭПЛ, рассчитанной полуэмпирическими методами, соответствующие модели взаимодействия АММ которой, как правило, описываются в рамках некулоновской ЭЭВ с положением центров зарядов в точках, не совпадающих с ядрами. Но эта особенность полуэмпирической ЭПЛ не является принципиальной с точки зрения точности, что было показано в работе [100] с использованием системы “атомных” центров, отличных от ядер, например, областей, отвечающих обращению в ноль их диполей.

Предпочтительность возможностей IEHT относительно CNDO-схемы на уровне описания атомных квадруполей [99] может представлять интерес для современных работ, посвященных мультипольному представлению полуэмпирической ЭПЛ с целью молекуляно-механического описания сложных систем в рамках КМ/ММ-подходов. Заслуживает упоминания использование приближения NDDO с базисами, ограниченными s-, p- и d-функциями и оценками АММ до квадруполя включительно [101]. Показано, что полученные ЭПО коррелируют с ЭПО, рассчитанным МР2- (метод Меллера–Плессе второго порядка) и B3LYP-методами независимо от выбора способа расчета ЭПО по тестовому набору точек на поверхности Конолли.

Такой же верхний порядок АММ до квадруполя оценен в работе, развивающей собственное представление геминалей, построенных из АО функций s- и p-типов также с использованием серии полуэмпирических (AM1, PM3, MNDO, NDDO) приближений [102, 103]. Для неполярных молекул СnН2n+ 2 метод с применением АММ в NDDO-версии показал себя как линейно масштабируемый, т.е., характеризуемый линейным ростом времени расчета при увеличении молекулы. Возможно, что необходимость точного учета дальнодействующих членов (на расстоянии более 20 Å) в более полярных молекулах такой показатель незначительно ухудшит, но данная версия обеспечивает отсутствие систематической ошибки в оценке теплот образования упомянутых молекул относительно ССП MNDO версии и выигрыш более чем на 2 порядка в скорости по сравнению с локальным ССП подходом. Последний метод в геминальном представлении был применен к описанию АММ одномерных цепочек (Н, С атомов и LiH) [104] и представляет интерес для полупроводниковых кристаллов.

1.2. Методы распределенного мультипольного анализа для твердотельных систем

1.2.1. Метод РМА по Саундерсу. Схема распределенного мультипольного анализа (РМА) Саундерса и колл. [1] является продолжением малликеновской схемы распределения ЭПЛ. Она отвечает разделению зарядовой плотности, образованной перекрыванием орбиталей соседних атомов, пополам между двумя соседними атомными центрами [36–39]. В рамках метода, используемого в пакете CRYSTAL и в данной работе, выражение для AMM $Q_{L}^{m}({\text{A}})$ данного кристаллографически независимого атома А в периодической бесконечной решетке имеет вид [1]:

(12)
$\begin{gathered} Q_{L}^{m}(A) = \sum\limits_\lambda {\sum\limits_1 {\sum\limits_2 {\sum\limits_g {P_{{12}}^{g}} } } } \times \\ \times \;\int {\chi _{1}^{0}(r)\chi _{2}^{g}(r)N_{L}^{m}X_{L}^{m}(r,A)dr,} \\ \end{gathered} $
где L и m являются, соответственно, порядком гармоники и проекции AMM, $P_{{12}}^{{\text{g}}}$ – элемент матрицы плотности, $\chi _{i}^{j}$ – i-атомная орбиталь в ячейке, смещенной на j-вектор, коэффициенты $N_{L}^{m} = (2 - {{\delta }_{{m0}}}) \times (L - \left| m \right|){\kern 1pt} !{\text{/}}(L + \left| m \right|){\kern 1pt} !$, ${{\delta }_{{m0}}}$ – дельта-функция, выражения ненормированных полиномов $X_{L}^{m}(r,Ai)$ Лежандра $X_{L}^{m}(r,A)$:
(13)
$\begin{gathered} X_{L}^{m}(r,A) = {{\Sigma }_{{(t,u,v)}}}{{D}_{L}}(t,u,v) \times \\ \times \;{{(x - {{X}_{{\text{A}}}})}^{t}}{{(y - {{Y}_{{\text{A}}}})}^{u}}{{(y - {{Z}_{{\text{A}}}})}^{v}} \\ \end{gathered} $
соответствуют ненормированным полиномам $X_{L}^{m}({\text{A}},i)$ Лежандра, взятых из [60] в декартовом представлении. Более сильная локализация АММ в твердом теле [40] существенно компенсирует тот недостаток, что в методе Саундерса и колл. затруднено применение слишком диффузных функций, которые приводят к катастрофически большому числу областей перекрывания и сбою. В силу использования АММ в пакете CRYSTAL как вспомогательного инструмента для замены двухэлектронных интегралов на одноэлектронные, в литературе редко рассматривалось поведение АММ относительно качества базиса [63, 105–107], типа DFT-функционала [63, 107]. Чаще изучалось влияние выбора базиса на заряды, заселенности связей [108]– [112] и ЭПО [105–107]). Оценивалось влияние на заряды поправок в дисперсионную ЭВ [112], типа зарядового распределения [112]. Более подробно изучалась геометрия и относительная стабильность полиморфных форм цеолитов в расчете на единицу SiO2 [63, 105–110, 112, 113].

1.2.2. Метод КАММ. Известно применение метода КАММ к цеолитам для выяснения роли электростатики в понижении барьеров активации процессов протонирования пропилена [114], ацетилена (с образованием –О–СН=СН2) [115] при замещении Si/Al на уровне кластерной модели размеров 2Т (при переходе от активного центра H3Si–O–SiH2–OH к H3Si–OH–AlH2–OH) [115] или 3Т (при переходе от активного центра Х3Si–O–SiH2–O–SiХ3 к H3Si–OH–AlH2–O–SiХ3, Х = = Н, ОН) [114], вписанным в разных положениях в решетке ZSM-5. Были учтены АММ до квадруполей выделенной зоны цеолита [115], взаимодействующей с активным центром, либо до суммарного четвертого порядка обоих атомов [114]. Однако, размеры зоны не уточнялись. Без учета вклада периодического кристаллического поля в точках с максимальным понижением барьеров активации (на стенке основного канала) электростатические взаимодействия доминировали в величине уменьшения барьеров.

1.2.3. Метод Хансена–Коппенса. Извлечение АММ из данных РС анализа по формализму Хансена–Коппенса с помощью пакетов MOLLY [116], XD [117], было продемонстрировано для цеолита натролита [118], берлинита [119], алюмодигидрофосфата аммония [120], Ala–Ser–Ala-трипептидов [121]. Предшественником MOLLY можно считать пакет VALRAY Стюарта [122], базирующийся на собственном построении АММ автором [123]. Авторы работы [121] пришли к выводу о возможной переносимости АММ (до октуполя) с лучшими перспективами для углеродных атомов основной цепи, но менее надежной для N- и Н-атомов, из которых проблема оценки Н-параметров по деформационной плотности РСА традиционна, а, кроме того, оба N- и Н-атома участвуют в ВС. Заметим, что АММ водородов Н(С) оценивались до диполей, кроме Н(N) амино-группы, где рассматривались и квадрупольные компоненты. Также была показана зависимость зарядов (С, N, О, Р) и компонентов октуполя атома С в серине от изменения геометрии [121]. Изменение атомных зарядов и доминирующей компоненты октупольного момента атома α‑углерода серина в пептидной цепи трех фрагментов аланин-серин-аланин было проанализировано относительно торсионного угла по оси α‑С–С', находящегося между α-С и атомом азота в серине [121]. Изменение зарядов α-С и С' указывает на важный перенос заряда, вариация компонент октуполя остается между 12 и 7% в интервале угла от –50 до 150°.

Сложности сравнения АММ, оцененных по Хансену–Коппенсу и АММ в рамках других схем, вытекают из зависимости коэффициентов сжатия-расширения сферической (κ) и асимметричной (κ') компонент ЭПЛ валентных орбиталей отдельных атомов от зарядного состояния атома. Общие трудности их оценки (κ, κ') обсуждались в работе [124]. По сравнению с АММ, рассчитанными из традиционных слейтеровских АО, изменение квантового числа nl (показатель радиальной переменной для слейтеровской АО) в обработке [116] описывается эффективными величинами для элементов разных периодов, которые приводят к дополнительному масштабированию подгоняемых АММ. В применении метода Хансена–Коппенса рассматриваются, как правило, наиболее значимые компоненты. Например, для описания ЭПЛ в трипептидах на С(sp2) имеют значение лишь компоненты P20, P33+, а для С(sp3) – P31–, P31+, P33+ [125]. Как результат, величины АММ имеют эффективное значение, и физический смысл имеет величина ЭПЛ, соответствующая данному набору АММ. Поэтому сравнение рассчитанных АММ с экспериментальными, полученными из подгонки данных РСА [116], возможно в смысле относительных величин АММ после перехода к атомной системе координат, принятых в РСА и отличающейся от таковой в РМА.

При сравнении экспериментальных атомных зарядов, объемов и дипольного момента п-нитроанилина, полученные с пакетами TOPXD [124], с рассчитанными на уровне ХФ и DFT, была показана необходимость κ-масштабирования теоретических величин и для ХФ, и для DFT-величин для того, чтобы достичь согласия.

1.2.4. Fast multipole method (FMM) или быстрый мультипольный метод. Расчеты разными вариантами метода fast multipole method (FMM) [126–129] ориентированы на полимеры [126], нанотрубки и кристаллы минералов [130], дендримеры [131]. Этот метод можно пояснить, рассматривая упрощенное представление системы в виде нескольких зон ячеек, в каждой из которых рассчитываются мультиполи ячеек. Размер ячеек каждой зоны определяется расстоянием от ее центра до места расчета ЭПО. Ячейки ближайшей зоны к месту расчета ЭП необходимо уменьшать разбиением на достаточно мелкие, взаимодействие с которыми учитывается более точно. Удаленные ячейки выбираются более крупными по оценке достаточной скорости сходимости. Разработаны разные варианты FMM, такие как ячеечный (cell FF или CFF) [126, 131, 132], который был ориентирован на распределения мультиполей (до квадруполей), возникающих при эмпирическом моделировании систем, и методы для мультиполей, появляющихся в DFT-расчетах [3, 130]. Благодаря линейной масштабируемости расчетов (рост времени счета линеен относительно числа атомов ∼ N), FMM-методы постепенно заменяют метод Эвальда (рост времени ∼ NlnN). В качестве примера в периодическом варианте pGvFMM [130] рассматривалась NaCl c суперъячейками до 512 атомов и нанотрубки (10, 10). На сегодняшний день нам известен единственный случай применения FMM к цеолитам [133].

2. ПРИМЕНЕНИЕ РАСПРЕДЕЛЕННОГО МУЛЬТИПОЛЬНОГО АНАЛИЗА К ОПИСАНИЮ ЭПО И ЭП В МОЛЕКУЛЯРНЫХ СИТАХ. ВОЗМОЖНОСТЬ ПОДСТРОЙКИ ЭПО ПОД НОВЫЙ БАЗИС

Развиваемый подход к описанию электростатического потенциала (ЭПО), поля (ЭП) и более высоких производных ЭПО в произвольных цеолитных системах (часть 1.1.2, уравнения (6), (7)) [63] акцентирован на количественном или полуколичественном рассмотрении, что требует оценку погрешности при выбранном базисе относительно расчета более высокого уровня. Поэтому должно быть изучено поведение ЭПО и его производных (ЭП, градиента ЭП) относительно повышения уровня базиса при аккуратном периодическом расчете на уровне функционала плотности. Нами были рассчитаны величины ЭПО и ЭП при разных базисах возрастающего качества, пронумерованных от 1 до 14 (таблица 1). Сравнение выполнялось на единой сетке точек, независимой от базиса, в трех разных сечениях, в которых доля доступного для адсорбата пространства была больше 50%. Формально граница контролировалась по линии нулевого потенциала, которая в наших системах хорошо воспроизводит “границу” полости цеолита и меняется с каждым базисом или уровнем расчета (рис. 1). Изменения доли доступного пространства при смене базиса приводят к относительно небольшим изменениям по сравнению с изменениями, которые потребовались бы при работе с сетками, связанными со структурой, оптимизированной с каждым использованным базисом.

Рис. 1.

Изоконтурные линии ЭПО = 0, отвечающие расчету на B3LYP уровне в сечении по атомам O37–O41–O45 цеолита MgPHI с базисами BS2, BS9, BS11 и BS14 (табл. 1) показаны точка-точка-штриховой, точечно-штриховой, короткой штриховой и штриховой линиями соответственно.

Чтобы описать сходимость ЭП и ЭПО с качеством базиса сравнивали относительную среднеквадратичную разницу (ОМНК или RRMS (метод относительных наименьших квадратов)) между сечениями ЭПО или ЭП (абсолютными величинами), рассчитанными с текущим BSN- и наиболее полным BSL-базисами:

(14)
$RRMS(EP_{{BSN}}^{{BSL}}) = \sqrt {\frac{{\sum\limits_{i = 1}^M {{{{(EP_{i}^{{BSN}} - EP_{i}^{{BSL}})}}^{2}}} }}{{\sum\limits_{i = 1}^M {{{{(EP_{i}^{{BSL}})}}^{2}}} }}} ,$
где М – число точек доступного пространства из общего числа точек сетки М0 (М0 ~ 104).

На примере серии объектов (для ПКЦ [61, 106], Н-форм [105, 107], катионных форм [106] были обнаружены: 1) сходимость ЭП с качеством базиса (рис. 2), 2) линейная зависимость ионности q0 (средний заряд О атомов) и ЭПО (рис. 3). Более подробного рассмотрения заслуживают результаты для катионной формы MgPHI. Координаты каркаса были фиксированы по результатам оптимизации c силовым полем Катлоу и др. [134] и пакетом GULP [135] для точного совпадения сеток, на которых рассчитываются ЭПО и остальные характеристики. Данные для всех 14 базисов (за исключением трех первых) описываются линейными зависимостями с коэффициентом корреляции более 0.9. На рис. 3 приведены данные для ПКЦ PHI (незаштрихованные символы), у которых наблюдается похожее поведение в данных координатах. Косвенным подтверждением удачного выбора сечения ЭПО = 0 оказалось то, что обработка данных в координатах q0(EP) (электростатический потенциал) приводит к линейным зависимостям с одинаковым наклоном в разных проекциях одного цеолита [105, 107]. Наклоны q0(EP) для разных цеолитов визуально отличаются.

Рис. 2.

Относительная среднеквадратичная разница (в долях) ЭП как функция типа базиса (табл. 1) в трех сечениях MgPHI (1–3 слева направо – кружки, треугольники, ромбы), показанных на вставке.

Рис. 3.

Зависимость ионности q0 (е) и ЭПО (a.e.) на B3LYP уровне с базисами, показанными в табл. 1, в трех сечениях, показанных на рис. 2 (данные в сечениях 1, 2, 3 даны кружками, треугольниками, ромбами) для цеолита филлипсита (PHI) MgPHI (заштрихованные) и ПКЦ (полностью кремниевый цеолит) PHI (незаштрихованные). Сплошной, штриховой и точечной линиями показаны аппроксимации для ЭПО в MgPHI. Данные, выпадающие из линейного поведения, обведены прямоугольниками с указанием базисов.

В то же время сходимость ЭПО, в отличие от ЭП, с качеством базиса не достигается (рис. 4а). Сходимость ЭП (рис. 2) означает, что кривизна ЭПО поверхности сходится относительно качества базиса. Тогда при смене базиса новая ЭПО поверхность перемещается на некоторую почти постоянную величину, обусловленную скачком ионности. Если это так, то близость ЭПО поверхностей может быть проверена линейным преобразованием, которое вытекает из зависимости ионности и ЭПО (рис. 3). Мы оценили точность такого перемасштабирования ЭПО на примере сечений MgPHI по следующей схеме. Величины $EP_{{i,cal}}^{{BSL}}$, рассчитанные с базисом BSL, в каждой точке (xi, yi) были пересчитаны в соответствии со сдвигом от величины ионности $q_{0}^{{BSL}}$ к $q_{0}^{{BSN}}$, отвечающей базису BSN по выражению, следующему из линейной зависимости ЭП от ионности:

(15)
$EP_{{i,app}}^{{BSN}}({{x}_{i}},{{y}_{i}}) = EP_{{i,cal}}^{{BSL}}({{x}_{i}},{{y}_{i}}) + \frac{1}{a}(q_{0}^{{BSN}} - q_{0}^{{BSL}}).$
Рис. 4.

Относительная среднеквадратичная разница ЭПО (RRMS(EP), %) ЭПО (14) как функция типа базиса (табл. 1) до преобразования (15) (а) и после (б).

По такому принципу получены величины в строке “EР (4)” в табл. 7 работы [106]. Из сравнения рассчитанного ЭПО (EР(2)) и скорректированного ЭПО (EР (4)) в табл. 7 работы [106] можно обнаружить близость соответствующих ОМНК, начиная с базиса BS4, с которого все базисы подчиняются линейной зависимости “ионность – ЭПО” (рис. 4б). Как правило, чем больше разница ионностей между двумя базисами, тем больше уменьшение ошибки (14) за счет сдвига ЭПО (15). Например, разности ионностей максимальны между BS14 и BS9 (0.507 e) или BS14 между и BS7 (0.282 e). Соответствующее уменьшение ОМНК (или RRMS) ЭПО за счет сдвига (15) может быть примерно оценено как от 27 до 3% для BS7 и от 47–50% до 6% для BS9. Для остальных базисов ОМНК ЭПО уменьшается до величины, сопоставимой с ОМНК ЭП на уровне расчета с тем же базисом.

Такое преобразование позволяет оценить в ОМНК разнице двух ЭПО долю, обусловленную различием ионностей, рассчитываемых с каждым из базисов. Оставшаяся ошибка определяется, в основном, различием кривизны ЭПО, т.е., различием ЭП. Это подтверждается сравнением ОМНК ЭПО верхней (EF (электрическое поле) (2)) и нижней (EP(4)) строк табл. 7 работы [106]. Для всех базисов после BS4 нaблюдаются близкие величины ОМНК ЭП и ЭПО. Дополнительным фактором различий между сдвинутой ЭПО поверхностью и наиболее точной является перераспределение зарядов и АММ при изменении ионности. Такая задача может быть решена с применением схемы (6), (7) . Однако, поскольку при всех BSN-базисах c N > 4 разница сохранялась менее 10%, то мы сочли такое совпадение удовлетворительным и не добивались более точного.

ВЫВОДЫ

Устранить переопределенность АММ в рамках рассмотренных выше методов можно было бы с использованием корреляций между зарядами и АММ соседних атомов в форме (6), (7). Для описания АММ произвольного цеолита предложена компактная форма из 8N' параметров (при учете до гексадекапольных моментов включительно), где N' – число типов атомов (Si, Al, O, …), из которых можно определить все m-компоненты АММ порядка L. Число 8N' определяется из двух параметров ${{a}_{L}}~\;{\text{и}}\;{\text{ }}{{b}_{L}}$ (6) на каждый АММ до порядка L = 4. Еще 6N' потребуется на описание всех зарядов [61, 62]. Данный набор заменяет $N\sum\nolimits_{L = 1}^4 {(2L + 1)} $ величин АММ (и еще N для зарядов), где (2L + 1) проекций моментов для каждого из L моментов всех $N = \sum\nolimits_{j = 1}^{N'} {{{N}_{j}}} $ атомов цеолита в ячейке, Nj – число кристаллографически независимых атомов каждого типа. Более того, каждый из 14N' параметров, как правило, применим и для оценки АММ таких же атомов другой системы того же или близкого состава, включая их искаженные каркасы, например, в случае колебаний. Любая комбинация компонент АММ L-порядка может быть аппроксимирована общей функцией вида (6), даже если точные величины для такой комбинации не рассчитываются пакетом CRYSTAL. Эта позволяет рассчитать все компоненты квадрупольного и октупольного моментов в том виде, в котором они используются в пакете GAMESS вместо того представления, в котором эти компоненты можно получить с помощью CRYSTAL. Решение данной проблемы требуется, например, для оценки всех десяти компонент октупольного момента, используемых в пакете GAMESS, вместо семи, которые можно получить с помощью CRYSTAL и аппроксимировать, как описано в [61–66, 105–107].

Результаты процедуры сдвига ЭПО (15) означают, что отсутствие сходимости ЭПО относительно улучшения базиса в ситах можно исправить с учетом изменения ионности, так что достигается минимальное отклонение ЭПО от наиболее точного, лимитируемое различием ЭП. Данное подобие позволяет преобразовывать трехмерную сетку сечений ЭПО к ионности, отвечающей нужному базису, например, который применяется для моделирования центрального кластера на КМ-уровне. Для этого достаточно оценить ионность, соответствующую нужному базису. Последняя задача может быть решена с помощью аппроксимации зарядов части/всех элементов цеолита в рамках выражений, полученных в [61–64]. После этого можно использовать выражение (15). При этом сходимость векторного потенциала электрического поля в положении центрального кластера на КМ-уровне обеспечивается выбором базиса для окружающей оболочки вокруг центрального кластера из широкого набора базисов, которые характеризуются сходимостью. В нашем случае (табл. 1), это базисы, более высокого качества, чем BS3, для которых ОМНК(ЭП) менее 10% (рис. 2).

Список литературы

  1. Saunders V.R., Freyria-Fava C., Dovesi R. et al. // Mol. Phys. 1992. V. 77. № 4. P. 629.

  2. Frisch M.J.G., Trucks W., Schlegel H.B. et al. // Gaussian 16. Inc.: Wallingford, CT. 2016.

  3. Strain M.C., Scuseria G.E., Frisch M.J. // Science. 1996. V. 271. № 5245. P. 51.

  4. Dovesi R., Orlando R., Erba A. et al. // Int. J. Quantum Chem. 2014. V. 114. № 19. P. 1287.

  5. Tomasi J., Mennucci B., Cammi R. // Chem. Rev. 2005. V. 105. № 8. P. 2999.

  6. Stone A.J. // Chem. Phys. Lett. 1981. V. 83. № 2. P. 233.

  7. Freitag M.A., Gordon M.S., Jensen J.H. et al. // J. Chem. Phys. 2000. V. 112. № 17. P. 7300.

  8. Granovsky A.A. // www http//classic.chem.msu.su/ gran/firefly/index.html. 2014.

  9. Schmidt M.W., Baldridge K.K., Boatz J.A. et al. // J. Comput. Chem. 1993. V. 14. № 11. P. 1347.

  10. Hall G.G. // Chem. Phys. Lett. 1973. V. 20. № 6. P. 501.

  11. Hall G.G. // Int. J. Quantum Chem. 1975. V. 9. № S9. P. 279.

  12. Martin D., Hall G.G. // Theor. Chim. Acta. 1981. V. 59. № 3. P. 281.

  13. Bonaccorsi R., Scrocco E., Tomasi J. // J. Am. Chem. Soc. 1976. V. 98. № 14. P. 4049.

  14. Bonaccorsi R., Scrocco E., Tomasi J. // Ibid. 1977. V. 99. № 14. P. 4546.

  15. Brobjer J.T., Murrell J.N. // J. Chem. Soc. Faraday Trans. 2 Mol. Chem. Phys. 1982. V. 78. № 11. P. 1853–1870.

  16. Bentley J. // Chemical Applications of Atomic and Molecular Electrostatic Potentials: Reactivity, Structure, Scattering, and Energetics of Organic, Inorganic, and Biological Systems / ed. Politzer P., Truhlar D.G. Boston, MA: Springer US, 1981. P. 63.

  17. Shipman L.L. // Chem. Phys. Lett. 1975. V. 31. № 2. P. 361.

  18. Tough R.J.A., Stone A.J. // J. Phys. A. Math. Gen. 1977. V. 10. № 8. P. 1261.

  19. Whitten A.E., Radford C.J., McKinnon J.J. et al. // J. Chem. Phys. 2006. V. 124. № 7. P. 074106.

  20. Bultinck P., Van Alsenoy C., Ayers P.W. et al. // Ibid. 2007. V. 126. № 14. P. 144111.

  21. Lillestolen T.C., Wheatley R.J. // Chem. Commun. 2008. V. 7345. № 45. P. 5909.

  22. Lillestolen T.C., Wheatley R.J. // J. Chem. Phys. 2009. V. 131. № 14.

  23. Manz T.A., Sholl D.S. // J. Chem. Theory Comput. 2010. V. 6. № 8. P. 2455.

  24. Manz T.A., Sholl D.S. // Ibid. 2012. V. 8. № 8. P. 2844.

  25. Geldof D., Krishtal A., Blockhuys F. et al. // J. Chem. Phys. 2014. V. 140. № 14.

  26. Geldof D., Krishtal A., Geerlings P. et al. // J. Phys. Chem. A. 2011. V. 115. № 45. P. 13096.

  27. Guerra C.F., Handgraaf J.-W., Baerends E.J. et al. // J. Comput. Chem. 2004. V. 25. № Vdd. P. 189.

  28. Rousseau B., Peeters A., Van Alsenoy C. // J. Mol. Struct. THEOCHEM. 2001. V. 538. № 1–3. P. 235.

  29. Rafat M., Popelier P.L.A. // J. Chem. Phys. 2006. V. 124. № 14. P. 144102.

  30. Joubert L., Popelier P.L.A. // Mol. Phys. 2002. V. 100. № 21. P. 3357.

  31. Popelier P.L.A., Rafat M. // Chem. Phys. Lett. 2003. V. 376. № 1–2. P. 148.

  32. Laidig K.E. // J. Phys. Chem. 1993. V. 97. № 49. P. 12760.

  33. Bader R.F.W., Keith T.A., Gough K.M. et al. // Mol. Phys. 1992. V. 75. № 5. P. 1167.

  34. Kosov D.S., Popelier P.L.A. // J. Phys. Chem. A. 2000. V. 104. № 31. P. 7339.

  35. Löwdin P. // J. Chem. Phys. 1950. V. 18. № 3. P. 365.

  36. Mulliken R.S. // Ibid. 1955. V. 23. № 10. P. 1833.

  37. Mulliken R.S. // Ibid. 1955. V. 23. № 10. P. 1841.

  38. Mulliken R.S. // Ibid. 1955. V. 23. № 12. P. 2338.

  39. Mulliken R.S. // Ibid. 1955. V. 23. № 12. P. 2343.

  40. Stone A. The Theory of Intermolecular Forces. Oxford University Press, 2013.

  41. Reed A.E., Weinstock R.B., Weinhold F. // J. Chem. Phys. 1985. V. 83. № 2. P. 735.

  42. Reed A.E., Curtiss L.A., Weinhold F. // Chem. Rev. 1988. V. 88. № 6. P. 899.

  43. Popelier P.L.A. // Phys. Scr. IOP Publishing, 2016. V. 91. № 3. P. 033007.

  44. Spackman M.A. // J. Chem. Phys. 1986. V. 85. № 11. P. 6587.

  45. Sokalski W.A., Sawaryn A. // J. Chem. Phys. 1987. V. 87. № 1. P. 526.

  46. Cooper D.L., Stutchbury N.C.J. // Chem. Phys. Lett. 1985. V. 120. № 2. P. 167.

  47. Whitehead C.E., Breneman C.M., Sukumar N. et al. // J. Comput. Chem. 2003. V. 24. № 4. P. 512.

  48. Popelier P. // Theor. Chim. Acta. 2001. V. 105. P. 393.

  49. Stone A.J. // J. Chem. Theory Comput. American Chemical Society, 2005. V. 1. № 6. P. 1128.

  50. Verstraelen T., Ayers P.W., Van Speybroeck V. et al. // Chem. Phys. Lett. Elsevier B.V., 2012. V. 545. P. 138.

  51. Misquitta A.J., Stone A.J., Fazeli F. // J. Chem. Theory Comput. 2014. V. 10. № 12. P. 5405.

  52. Koch U., Egert E. // J. Comput. Chem. 1995. V. 16. № 8. P. 937.

  53. Price S.L., Stone A.J. // J. Chem. Soc. Faraday Trans. 1992. V. 88. № 13. P. 1755.

  54. Stone A.J., Alderton M. // Mol. Phys. 2002. V. 100. № 1. P. 221.

  55. Sokalski W.A., Poirier R.A. // Chem. Phys. Lett. 1983. V. 98. № 1. P. 86.

  56. Stone A.J., Price S.L. // J. Phys. Chem. 1988. V. 92. № 12. P. 3325.

  57. Price S.L. // J. Chem. Soc. Faraday Trans. 1996. V. 92. № 17. P. 2997.

  58. Price S.L., Hamad S., Torrisi A. et al. // Mol. Simul. 2006. V. 32. № 12–13. P. 985.

  59. Wheatley R.J. // Chem. Phys. Lett. 1993. V. 208. № 3–4. P. 159.

  60. Pisani C., Dovesi R., Roetti C. Berlin, Heidelberg: Springer Berlin Heidelberg, 1988. V. 48.

  61. Larin A.V. // J. Comput. Chem. 2011. V. 32. № 11. P. 2459.

  62. Larin A.V. // J. Struct. Chem. 2014. V. 55. № 3. P. 398.

  63. Larin A.V., Parbuzin V.S., Vercauteren D.P. // Int. J. Quantum Chem. 2005. V. 101. № 6. P. 807.

  64. Larin A.V., Trubnikov D.N., Vercauteren D.P. // Russ. J. Phys. Chem. A. 2005. V. 79. № 6. P. 986.

  65. Larin A.V., Hansenne C., Trubnikov D.N. et al. // Int. J. Quantum Chem. 2005. V. 105. № 6. P. 839.

  66. Larin A.V., Semenyuk R.A., Trubnikov D.N. et al. // Russ. J. Phys. Chem. A. 2007. V. 81. № 4. P. 493.

  67. Becke A.D. // J. Chem. Phys. 1988. V. 88. № 4. P. 2547.

  68. Murray C.W., Handy N.C., Laming G.J. // Mol. Phys. 1993. V. 78. № 4. P. 997.

  69. Price S.L., Faerman C.H., Murray C.W. // J. Comput. Chem. 1991. V. 12. № 10. P. 1187.

  70. Faerman C.H., Price S.L. // J. Am. Chem. Soc. 1990. V. 112. № 12. P. 4915.

  71. Zakrzewska K., Pullman A. // J. Comput. Chem. 1985. V. 6. № 4. P. 265.

  72. Bellido M.N., Rullmann J.A.C. // J. Comput. Chem. 1989. V. 10. № 4. P. 479.

  73. Kędzierski P., Sokalski W.A. // J. Comput. Chem. 2001. V. 22. № 10. P. 1082.

  74. Sokalski W.A., Shibata M., Rein R. et al. // J. Comput. Chem. 1992. V. 13. № 7. P. 883.

  75. Maroulis G. // J. Phys. Chem. 1996. V. 100. № 32. P. 13466.

  76. Sokalski W.A., Maruszewski K., Hariharan P.C. et al. // Int. J. Quantum Chem. 1989. V. 36. № 16 S. P. 119.

  77. Jakobsen S., Kristensen K., Jensen F. // J. Chem. Theory Comput. 2013. V. 9. № 9. P. 3978.

  78. Jakobsen S., Jensen F. // J. Chem. Theory Comput. 2014. V. 10. № 12. P. 5493.

  79. Breneman C.M., Wiberg K.B. // J. Comput. Chem. 1990. V. 11. № 3. P. 361.

  80. Chirlian L.E., Francl M.M. // J. Comput. Chem. 1987. V. 8. № 6. P. 894.

  81. Singh U.C., Kollman P.A. // Ibid. 1984. V. 5. № 2. P. 129.

  82. Swart M., van Duijnen P.T., Snijders J.G. // J. Comput. Chem. 2001. V. 22. № 1. P. 79.

  83. Simmonett A.C., Gilbert A.T.B., Gill P.M.W. // Mol. Phys. 2005. V. 103. № 20. P. 2789.

  84. Devereux M., Raghunathan S., Fedorov D.G. et al. // J. Chem. Theory Comput. 2014. V. 10. № 10. P. 4229.

  85. Momany F.A. // J. Phys. Chem. 1978. V. 82. № 5. P. 592.

  86. Bayly C.I., Cieplak P., Cornell W. et al. // J. Phys. Chem. 1993. V. 97. № 40. P. 10269.

  87. Campañá C., Mussard B., Woo T.K. // J. Chem. Theory Comput. 2009. V. 5. № 10. P. 2866.

  88. Wiberg K.B., Rablen P.R. // J. Comput. Chem. 1993. V. 14. № 12. P. 1504.

  89. Cieplak P., Dupradeau F.Y., Duan Y. et al. // J. Phys. Condens. Matter. 2009. V. 21. № 33.

  90. Tsiper E.V., Burke K. // J. Chem. Phys. 2004. V. 120. № 3. P. 1153.

  91. Williams D.E. // J. Comput. Chem. 1988. T. 9. № 7. C. 745.

  92. Nistor R.A., Polihronov J.G., Müser M.H. et al. // J. Chem. Phys. 2006. V. 125. № 9.

  93. Verstraelen T., Sukhomlinov S.V., Van Speybroeck V. et al. // J. Phys. Chem. C. American Chemical Society, 2012. V. 116. № 1. P. 490.

  94. Hirshfeld F.L. // Theor. Chim. Acta. 1977. V. 44. № 2. P. 129.

  95. Parr R.G., Ayers P.W., Nalewajski R.F. // J. Phys. Chem. A. 2005. V. 109. № 17. P. 3957.

  96. Bertoni C., Slipchenko L.V., Misquitta A.J. et al. // Ibid. 2017. V. 121. № 9. P. 2056.

  97. Verstraelen T., Vandenbrande S., Heidar-Zadeh F. et al. // J. Chem. Theory Comput. 2016. V. 12. № 8. P. 3894.

  98. Misquitta A.J., Stone A.J. // Theor. Chem. Acc. Springer Berlin Heidelberg, 2018. V. 137. № 11. P. 153.

  99. Rein R. // Advances in Quantum Chemistry. 1973. V. 7. № C. P. 335.

  100. Dovesi R., Pisani C., Ricca F. et al. // J. Chem. Soc. Faraday Trans. 2. 1974. V. 70. P. 1381.

  101. Horn A.H.C., Lin H., Clark T. // Theor. Chem. Acc. 2005. V. 114. № 1–3. P. 159.

  102. Tokmachev A.M., Tchougréeff A.L. // J. Phys. Chem. A. 2005. V. 109.№ 33. P. 7613.

  103. Tokmachev A.M., Tchougréeff A.L. // Ibid. 2003. V. 107. № 3. P. 358.

  104. Tokmachev A.M., Dronskowski R. // Chem. Phys. 2006. V. 322. № 3. P. 423.

  105. Larin A.V., Trubnikov D.N., Vercauteren D.P. // J. Comput. Chem. 2008. V. 29. № 1. P. 130.

  106. Larin A.V., Sakodynskaya I.K., Trubnikov D.N. // Ibid. 2008. V. 29. № 14. P. 2344.

  107. Larin A.V., Trubnikov D.N., Vercauteren D.P. // Int. J. Quantum Chem. 2007. V. 107. № 15. P. 3137.

  108. Aprà E., Dovesi R., Freyria-Fava C. et al. // Model. Simul. Mater. Sci. Eng. 1993. V. 1. № 3. P. 297.

  109. White J.C., Nicholas J.B., Hess A.C. // J. Phys. Chem. B. 1997. V. 101. № 4. P. 590.

  110. White J.C., Hess A.C. // J. Phys. Chem. 1993. V. 97. № 24. P. 6398.

  111. Zicovich-Wilson C.M., Hô M., Navarrete-López A.M. et al. // Theor. Chem. Acc. 2016. V. 135. № 8. P. 188.

  112. Román-Román E.I., Zicovich-Wilson C.M. // Chem. Phys. Lett. Elsevier B.V. 2015. V. 619. P. 109.

  113. Civalleri B., Zicovich-Wilson C., Ugliengo P. et al. // Chem. Phys. Lett. 1998. V. 292, № 4–6. P. 394.

  114. Dziekoński P., Sokalski W.A., Szyja B. et al. // Ibid. 2002. V. 364. № 1–2. P. 133.

  115. Dziekoński P., Sokalski W.A., Kassab E. et al. // Ibid. 1998. V. 288. № 2–4. P. 538.

  116. Hansen N.K., Coppens P. // Acta Crystallogr. Sect. A. 1978. V. 34. № 6. P. 909.

  117. Volkov A., Macchi P., Farrugia L.J., et al. XD2006 − A Computer Program Package for Multipole Refinement, Topological Analysis of Charge Densities and Evaluation of Intermolecular Energies from Experimental and Theoretical Structure Factors. 2006.

  118. Ghermani N.E., Lecomte C., Dusausoy Y. // Phys. Rev. B. 1996. V. 53. № 9. P. 5231.

  119. Aubert E., Porcher F., Souhassou M. et al. // Acta Crystallogr. Sect. B Struct. Sci. 2003. V. 59. № 6. P. 687.

  120. Pérès N., Boukhris A., Souhassou M. et al. // Acta Crystallogr. Sect. A Found. Crystallogr. 1999. V. 55. № 6. P. 1038.

  121. Koritsanszky T., Volkov A., Coppens P. // Ibid. 2002. V. 58. № 5. P. 464.

  122. Stewart R.F. // Acta Crystallogr. Sect. A. 1976. V. 32. № 4. P. 565.

  123. Stewart R.F. // Isr. J. Chem. 1977. V. 16. № 2–3. P. 124.

  124. Volkov A., Gatti C., Abramov Y. et al. // Acta Crystallogr. Sect. A Found. Crystallogr. 2000. V. 56. № 3. P. 252.

  125. Pichon-Pesme V., Lecomte C., Lachekar H. // J. Phys. Chem. 1995. V. 99. № 16. P. 6242.

  126. Ding H.Q., Karasawa N., Goddard W.A. // J. Chem. Phys. 1992. V. 97. № 6. P. 4309.

  127. Nijboer B.R.A., De Wette F.W. // Physica. 1957. V. 23. № 1–5. P. 309.

  128. Lambert C.G., Darden T.A., Board J.A. // J. Comput. Phys. 1996. V. 126. № 2. P. 274.

  129. Challacombe M., White C., Head-Gordon M. // J. Chem. Phys. 1997. V. 107. № 23. P. 10131.

  130. Kudin K.N., Scuseria G.E. // Chem. Phys. Lett. 1998. V. 289. № 5–6. P. 611.

  131. Zheng J., Balasundaram R., Gehrke S.H. et al. // J. Chem. Phys. 2003. V. 118. № 12. P. 5347.

  132. Ding H.Q., Karasawa N., Goddard W.A. // Chem. Phys. Lett. 1992. V. 196. № 1–2. P. 6.

  133. Jousse F., Auerbach S.M. // J. Chem. Phys. 1997. V. 107/ № 22. P. 9629.

  134. Schröder K.-P., Sauer J., Leslie M. et al. // Chem. Phys. Lett. 1992. V. 188. № 3. P. 320.

  135. Gale J.D. // J. Chem. Soc. Faraday Trans. The Royal Society of Chemistry. 1997. V. 93. № 4. P. 629.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал физической химии

Архивы выпусков Информация о журнале Отправить рукопись в журнал