CC BY
11
УДК 539.192
МАСШТАБИРОВАНИЕ КВАНТОВОМЕХАНИЧЕСКОГО СИЛОВОГО ПОЛЯ МОЛЕКУЛЫ ПРОТИВ РЕШЕНИЯ ОБРАТНОЙ КОЛЕБАТЕЛЬНОЙ ЗАДАЧИ Ю.Н. Панченко, Ж.Р. Де Марэ*
(кафедра физической химии)
Рассмотрены особенности, характеризующие традиционный подход в вычислительной колебательной спектроскопии (обратная колебательная задача) и подход, основанный на использовании масштабированных квантовомеханических силовых полей молекул. Обсуждены некоторые результаты определения равновесной геометрии бензола как в гармоническом приближении, так и в приближении, учитывающем кинематическую и динамическую поправки на ангармоничность посредством решения обратной колебательной задачи. Сопоставлены результаты определения масштабирующих множителей с помощью разной математической техники расчетов для квантовомеханических силовых полей молекулы С2Е6, рассчитанных в качестве примера в трех различных теоретических приближениях.
Колебательная спектроскопия представляет собой один из экспериментальных методов исследования строения молекулы. Изучение колебательных спектров молекул позволяет определить симметрию молекулы, а также установить присутствие определенных функциональных групп и решить ряд вопросов внутримолекулярной динамики. Вычислительная колебательная спектроскопия играет существенную роль при получении такого рода результатов. Успехи техники квантовомеханических расчетов [1] привели к отходу от традиционной идеологии вычислительной колебательной спектроскопии, основанной на решении обратной колебательной задачи [2-4]. Поэтому ниже рассмотрены основные черты как традиционной идеологии, так и новой идеологии, которая сводится к поправкам для квантовомеханических силовых полей с помощью эмпирических масштабирующих множителей.
Как традиционный, так и новый подходы основаны на получении рассчитанных частот колебаний и их колебательных мод (форм) с помощью векового уравнения теории малых колебаний [2-4], т.е.
Здесь С - матрица кинематических коэффициентов, Е - матрица силовых постоянных, Ь - матрица нормированных колебательных мод, а Л - диагональная матрица частотных параметров Хк. Они связаны с частотами гармонических колебаний ю к соотношением Хк = dюk, где коэффициент d зависит от выбора единиц. Все приведенные выше матрицы -квадратные порядка п, где п - число независимых внутренних координат (п = 3Ы-6 для нелинейных и 3N-5 для линейных молекул, где N - число атомов в молекуле). Обе матрицы С и Е симметричные, неособенные и положительно определенные. Колебательные моды нормированы таким образом, что выполняется следующее соотношение:
LFL = Л или Хк = ~k F!k.
(2)
GFL = LЛ.
(1)
Символ тильда (~) означает транспонирование. Здесь Хк и ! - частотный параметр и мода к-го нормального колебания соответственно (! - к-й столбец матрицы Ь).
Традиционный подход
При этом подходе элементы матрицы С можно определить однозначно для принятой структуры мо-
G. R. De Maré , Service de Chimie Physique et de Photophysique (Atoms, Molécules et Atmosphères), Faculté des Sciences, Université Libre de Bruxelles, CP 160/09, 50 avenue F.-D. Roosevelt, B-1050 Brussels, Belgium, e-mail address: [email protected]
лекулы. В этих расчетах обычно используют экспериментальные геометрические параметры (если они доступны). В традиционном подходе определение Ж основано на следующих двух предположениях.
1. Для двух данных колебательных координат соответствующая пара диагональных элементов матрицы Ж существенно больше, чем соответствующий недиагональный элемент, связанный с этими координатами.
2. Силовые постоянные переносимы в ряду родственных молекул. Небольшие различия между экспериментальными частотами колебаний и частотами, рассчитанными с матрицей Ж, составленной из силовых постоянных родственных молекул, можно устранить, решая так называемую обратную колебательную задачу. При этом исправленное силовое поле должно быть "близким" к исходному силовому полю.
Для малых молекул эти предположения были подтверждены практическими расчетами (начало 50-х -70-е гг. прошлого века) и, следовательно, они рассматриваются как аксиомы при традиционном подходе к вычислительной колебательной спектроскопии. Таким образом, традиционный подход к определению матрицы Ж в вычислительной колебательной спектроскопии сводится к решению обратной колебательной задачи [2-4].
Основные проблемы при решении этой задачи -отыскание экспериментальных данных, относящихся к физическим величинам, связанным с силовыми постоянными, и выбор одного из п! возможных отнесений экспериментальных частот колебаний. Дело в том, что число силовых постоянных, которые следует уточнить, обычно больше, чем число легко доступных экспериментальных данных. Вследствие этого и на основе первого предположения предложено большое число моделей силовых полей с нулевыми силовыми постоянными, соответствующими дальним "взаимодействиям". Например, ограничения, соответствующие наиболее известной модели модифицированного валентно-силового поля (МВСП), накладывались в [5] даже на квантовомеханическое силовое поле в приближении МР2/6-3Ш*.
Прежде всего вековое уравнение (1) решают с исходной матрицей Затем матрицу уточняют с помощью итерационной процедуры, при которой используют, как правило, только эксперимен-
тальные частотные параметры Лехр(1. Минимальная норма вектора отклонений (невязок) частотных параметров (рассчитываемая на каждом итерационном шаге) от экспериментальных значений, обозначенная как ДА, служит критерием оптимизации в простейшем случае. Пусть т будет числом частот, известных из эксперимента (т>п). Тогда размерность ДА равна т. Если, следуя обычной практике, используют эвклидову норму ||ДА|| = (дА да) 2 , то задача сводится к минимизации функционала нелинейного метода наименьших квадратов, т.е.
т , . -
Д(Я) = Д~ -ДЬ=£(-АТ )
(3)
Здесь 8 - вектор нескольких параметров, регулирующих значения искомых частотных параметров ДАк. В частности такими параметрами могут быть силовые постоянные 5, которые уточняются при итерационной процедуре (пусть 5 обозначает в этом разделе элемент матрицы Ж на каждом итерационном шаге), а А^Х - к-й частотный параметр, соответствующий к-му экспериментальному значению. Однако помимо учета близости рассчитанных и экспериментальных частот существуют другие способы параметризации задачи (1) и другие критерии для оптимального решения (см. ниже).
Эффективным методом минимизации функционалов типа (3), который обладает асимптотически квадратичной сходимостью, является итерационная процедура Гаусса-Ньютона [6, 7], основанная на локальной линеаризации ДА(8) в окрестности исходного приближения 80:
ДА(8) = ДА(80) + J Д8. (4)
Здесь J - матрица Якоби с элементами J к] = Э А/Эл ], рассчитанными для набора параметров (силовых постоянных) 80, а Д8 = 8 - 80 означает вектор поправок для перехода к следующему приближению, который находят, решая систему линейных уравнений:
JJ Д8 = ^ ДА(80).
(5)
Малость нормы ||Д8|| служит критерием для окончания итерационной процедуры.
Как известно, если обратная колебательная задача плохо обусловлена, то необходимо использовать метод регуляризации, т.е. некоторое дополнительное условие. С точки зрения традиционной идеологии чисто эмпирического подхода к решению обратной
колебательной задачи (см. предположение 2), логичным представляется введение специальной штрафной функции, которая дает решение, близкое к исходной матрице силовых постоянных. В качестве меры близости исходной и конечной матриц силовых постоянных в [8-10] выбрано минимальное значение эвклидовой нормы их разности. Таким образом, этот подход использует регуляризацию по Тихонову [9-11]. Действительно, степень отклонения конечной матрицы силовых постоянных Е от ЕщМа1 учитывается здесь явным образом с помощью следующего выражения:
М(8) = К(8) + а||Е-Е ||2.
(6)
Здесь Л(8) означает обычный функционал нелинейного метода наименьших квадратов, определенный в (3), а а > 0 - параметр регуляризации (см., например, [8-11]). Соотношение (6) определено как сумма квадратов норм в двух различных простран-
т
ствах: К , где т означает число экспериментальных частот, согласованное с принятой моделью силового поля, и К , где ^ = п(п + 1)/2 обозначает число силовых постоянных в Е или Едт (принимая во внимание симметрию Е и Е?_т). Тем самым квантово-механическое силовое поле Ед—т включено в решение обратной колебательной задачи как Е1пШа1. В последующем в качестве Е1п1(1а1 будем рассматривать только исходную квантовомеханическую матрицу силовых постоянных Е .
д—т
Относительные вклады каждого члена в уравнение (6) зависят от выбора а. Обычно используется довольно малое значение а [8-11], в частности из-за значительного различия между размерностями про-
т ?
странств К и К . Однако важно не только различие между пространствами К и К, но в большей степени также различие в природе X и Е и в различном порядке их величин. По существу, ищут компромисс между двумя требованиями (в общем несовместимыми): максимальная близость Е к Е согласно эвклидовой норме и наиболее близкое воспроизведение
„ л ехР
экспериментальных значений Хк при помощи метода наименьших квадратов. Этот выбор не противоречит идеологии построения метода наименьших квадратов для чисто эмпирического решения обратной колебательной задачи, где искомая матрица не должна существенно отличаться от исходной чисто эмпирической Е1п1(1а1, причем решение выбирается так,
чтобы АХ было достаточно малым. Этот алгоритм решения обратной колебательной задачи соответствует одной из хорошо известных процедур регуляризации для плохо обусловленных задач. Следует отметить, что использование функционала (6) для решения обратной колебательной задачи с исходным набором квантовомеханических силовых постоянных, исправленных различными весовыми множителями для каждой силовой постоянной, ранее предложено в [8]. Более того, штрафная функция, аналогичная включенной в функционал (6), широко применялась ранее также при эмпирической оптимизации геометрии молекул (см., например, [12]).
Численный пример решения обратной колебательной задачи
Техника нелинейного метода наименьших квадратов, обсужденная в [8-10] и основанная на идеологии решения обратной колебательной задачи [2-4], недавно использована для определения равновесной геометрии бензола [13]. Это сделано как в гармоническом приближении, так и в приближении, учитывающем кинематическую и динамическую ангармоничность бензола [13] с совместным использованием электронографических интенсивностей [14], частот колебаний легкой молекулы и ее изотопомера С6Б6 [15] и силового поля Ед—т, рассчитанного в приближении МР2/6-3Ш*//МР2/6-3Ш* [1]. В связи с кинематической ангармоничностью функционал (6) снабжен штрафной функцией в||Я - К1п1(1а1||2, где Я и Я;пМа1 - матрицы, связанные с геометрией молекулы, а в - постоянная, соответствующая относительным весам межъядерных расстояний, т.е.
М(8) = К(8) + а||Е-ЕдЛ2 + в|| Я - КШ1Л2. (6а)
Межъядерные равновесные расстояния г полученные в [13], представлены в табл. 1. Здесь приведены также равновесные расстояния для гармонического и ангармонического приближений двухатомных молекул СН (х2Пг) и С2 (х1^) [16]. Табл. 1 содержит также результаты определения геометрических параметров, полученных в [17] путем совместного анализа результатов квантовомеханического расчета на высоком уровне (метод связанных кластеров с учетом одно- и двукратных возбуждений и неитерационной поправки на трехкратные возбуждения в четырех экспонентном корреляционно-согласованном базисе, СС8Б(Т)/сс-рУ^2) и эмпирических значений
Сравнение равновесных расстояний ге(А) для формальных связей С-Н и С-С в молекулах СН, С2 и С6Н6
С-Н (х2Пг) СбНб [13] с,
Равновесное расстояние Гармоническое приближение [16] Ангармоническое приближение Гармоническое приближение Кинематическая ангармоничность Кинематическая + динамическая ангармоничность Гармс приб1
Формальная связь С-Н 1,1199 1,1304 1,0840 1,0914 1,0857 1,080:
с2 (x'Eg+) СбНб [13] сбнб
Гармоническое приближение [16] Ангармоническое приближение Гармоническое приближение Кинематическая ангармоничность Кинематическая + динамическая ангармоничность Гармс приб1
Формальная связь С-С 1,24253 1,24575 1,3951 1,3968 1,3929 1,391.
Рассчитано по данным в [16].
вращательных постоянных четырех изотопомеров бензола.
Данные табл. 1 показывают, что для двухатомных молекул СН (х2Пг) и С2 (х1^) удлинение связи вследствие ангармоничности предсказано равным около 0,01 и 0,003 А соответственно. Напротив, в [13] сообщено, что с применением регуляризирую-щего алгоритма, представленного функционалом (6) (см. также [8-10]), удлинение формальных связей С-Н в бензоле предсказано равным около 0,007 А (с учетом кинематической ангармоничности) и только около 0,002 А (с учетом кинематической и динамической ангармоничности). Поправки к межъядерным расстояниям С-Н в [13] для бензола, полученные с учетом только кинематической ангармоничности, приблизительно совпадают с соответствующей экспериментальной поправкой для межъядерного расстояния двухатомной молекулы СН, когда учитывают как кинематическую, так и динамическую ангармоничность [18] (табл. 1). Предсказанное уменьшение длины расстояния С-Н в бензоле [13] при переходе от структуры, учитывающей кинематическую ангармоничность, к структуре, принимающей во внимание как кинематическую, так и динамическую ангармоничность, вызывает сомнение. Поэтому поправки на кинематическую ангармоничность для расстояний С-Н, полученные в [13], представляются неправдоподобными. Действительно, для частот колебаний вклад кинематической ангармоничности обычно составляет приблизительно только 25% от полной поправки на ангармоничность (см., например, [18, 19]).
Длины формальных связей С-С для бензола, рассчитанные в [13], предсказывают удлинение за счет ангармоничности только в случае кинематической ангармоничности (около 0,002 А). Однако совместный учет кинематической и динамической ангармоничностей приводит к уменьшению параметра ге(С-С) приблизительно на 0,002 А [13]. Метод учета динамической ангармоничности, описанный в [20, 21], использован в [13] как для длин формальных связей С-Н, так и для связей С-С. Различные знаки в поправках на ангармоничность, полученные в [13] для этих межъядерных расстояний (см. табл. 1), указывают на некорректность расчетов, основанных на функционале (6а). Следовательно, результаты, полученные в [13], не соответствуют критерию удлине-
ния межъядерных расстояний вследствие учета ангармоничности. Необходимо отметить, что значения межъядерных расстояний г полученные в [13], отличны от результатов в [17]. Последние ниже приблизительно на 0,004 А как для С-Н, так и для С-С формальных связей. Это очевидно связано прежде всего с включением специальной штрафной функции Р||Я - КшШа1||2 в функционал (6а), использованный в [13]. Действительно квантовомеханичес-кое силовое поле, полученное в [13] в приближении МР2/6-3Ш*//МР2/6-3Ш*, приводит к отклонениям рассчитанных частот как для С-Н валентных колебаний (у20, е1и), так и для деформационных ко -лебаний кольца (у4, Ь2§) приблизительно на 200 см-1. Следовательно, требование близости искомого решения к согласно функционалу (6а) должно, несомненно, искажать конечное решение задачи и в частности колебательные моды.
В [13] результаты ге(С-С) = 1,397 А и ге(С-Н) = 1,087 А, полученные оптимизацией геометрии в приближении МР2/6-3Ю*, использованы в качестве КЫ(;а1 для штрафной функции. Напротив, в [17] использованы ге(С-С) = 1,3911 А и ге(С-Н) = 1,0800 А, полученные оптимизацией геометрии в приближении СС8Б(Т)/ссУ02. При этом в [17] подчеркнуто, что оптимизация геометрии на таком высоком уровне может дать теоретические межъядерные расстояния, точность которых лучше, чем 0,003 А. Это означает, что включение штрафной функции в функционал (6а) не может улучшить оптимизированные геометрические параметры, полученные в [13] на уровне МР2/6-3Ю*.
Важно также отметить, что в [13] частоты колебаний бензола и его изотопомера С6Б6 включены как дополнительные данные в задачу определения геометрических параметров по интенсивностям рассеяния [14]. Напротив, в [17] использованы эмпирические вращательные постоянные, которые более чувствительны к изменениям геометрических параметров молекулы, чем частоты колебаний. Заслуживает упоминания также то, что рассмотрение удлинения формальной связи вследствие ангармоничности согласуется с оценочными расстояниями ге в [17]. Следует отметить, что совместное варьирование геометрии и решение колебательной задачи в [13] означает выход из теоретического минимума на потенциальной поверхности, т.е. С и Е рассчитаны для различных геометрических параметров.
В общем, как отмечено ранее в [22], попытки вовлечь как можно больше дополнительных данных в такую обратную колебательную задачу могут дать отрицательный результат. Дело в том, что некоторые из этих дополнительных экспериментальных данных измерены с гораздо более низкой точностью, чем другие. Как правило, они основаны на различных физических моделях и определены как коэффициенты в разложениях, где членами более высокого порядка пренебрегают. Такое пренебрежение при определении параметров должно влиять на конечные рассчитанные значения, потому что эти параметры рассчитываются, а не измеряются непосредственно.
Идеология масштабирования квантовомеханического силового поля
Этот подход основан на оптимизации геометрии молекулы и квантовомеханическом расчете потенциальной энергии с последующим аналитическим определением и, наконец, масштабированием полученного квантовомеханического силового поля.
Основным предположением, использованным в этом подходе, является постоянство относительных ошибок при определении силовых постоянных для одинаковых структурных групп в ряду родственных молекул в том же самом квантовомеханическом приближении. Это предположение позволяет переносить наборы масштабирующих множителей, полученных для сравнительно малых хорошо изученных молекул, на значительно большие молекулы. Следовательно, постоянство относительных ошибок становится основным постулатом этого подхода.
Существуют различные подходы к масштабированию Е (см., например, [23]). Согласно современ-
q—m
ному методу Пулаи (см., например, [24-31]), Е во
q—m
внутренних координатах подвергается конгруэнтному преобразованию, т. е.
Е = ^1/2 Fq_mS1/2.
(7)
Здесь матрица Е рассчитана в приближении Хар-три-Фока, а 8 - диагональная матрица масштабирующих множителей. Некоторые преимущества метода масштабирования по Пулаи [24-31] по сравнению с другими методами масштабирования ранее рассмотрены в [22, 23, 28-32]. Однако наиболее важное преимущество метода Пулаи состоит в его теоретическом обосновании [27]. Очевидно, это мотивировало его широкое применение.
Метод Пулаи предназначен для введения эффективных эмпирических поправок на корреляционные эффекты в матрицу Е , полученную в хартри-фо-
q—m
ковском приближении. Систематическое завышение силовых постоянных, рассчитанных в хартри-фоков-ском приближении, требует таких поправок. Действительно, в [27] показано, что хартри-фоковские силовые постоянные приблизительно линейно связаны с точными силовыми постоянными, рассчитанными с тем же самым базисным набором в приближении конфигурационного взаимодействия. Это служит основанием для введения одного общего масштабирующего множителя. В случае нескольких масштабирующих множителей, что связано с использованием для их определения экспериментальных частот колебаний, недиагональные элементы матрицы Е
умножаются на ys,s,, т.е. связь становится не-
q—m к 1
линейной.
Матрица 8 имеет следующее строение:
S =Х slPl.
1=1
При этом считается, что п колебательных координат расклассифицированы на q групп. Каждая группа содержит эквивалентные или квазиэквивалентные (однотипные с точки зрения масштабирования) координаты. Р1 - диагональная матрица порядка п. Она содержит одинаковые значения, равные единице, в положениях, соответствующих эквивалентным, квазиэквивалентным или неэквивалентным координатам 1-й группы, и нули во всех других положениях. Обозначение s¡ теперь представляет собой масштабирующий множитель для этой группы эквивалентных, квазиэквивалентных или неэквивалентных координат.
Сумма всех матриц Р1 дает единичную матрицу, т.е.
q
Е р1=Е
1=1
В общем число масштабирующих множителей q может лежать в области между 1 и п. Однако учет эквивалентных и квазиэквивалентных координат обычно приводит к неравенству д<< п.
Если искомая матрица Е зависит от параметра s, тогда частотные параметры Лк и колебательные моды ! также зависят от этого параметра. Используя аналог теоремы Гельмана-Фейнмана [33, 34], производная частотного параметра Лк по s определяется по
соотношению (2), а именно
Э* = ~ (д¥/дя)£к
(8)
Соотношения (3), (4) и (5) также справедливы для определения масштабирующих множителей. Однако теперь 8 - вектор размерности д, который составлен из масштабирующих множителей */, участвующих в итерационной процедуре. Теперь в соотношениях (4) и (5) элементы матрицы Якоби J с элементами =
Э'Я /Э*. рассчитываются по соотношению (8), где
к j
производную легко получить дифференцированием уравнения (7), т.е.
№
2*1
2*1
-Р,
/ \1/2
= 2 Е - [-т Р + Р/ Г-т Р ].
2 £
Суммирование проводится по всем группам эквивалентных и квазиэквивалентных координат.
Следовательно, элементы принимают вид [35, 36]:
1 ~ I д
,
\1/2
Л,
Г Р + Р Г
, д-т / / д-т ,
уравнений (5) в произвольном случае, включающем присутствие вырождения, является разложение матрицы Якоби по сингулярным числам (см., например,
[40, 41]):
J = иX V, (10)
где матрицы и (порядка т) и V (порядка д) ортогональны, а матрица X (т строк, д столбцов) псевдо-диагональна с неотрицательными элементами X =о. Величины о - сингулярные числа J, а столбцы матрицы V - правые сингулярные векторы. Если между столбцами J существуют линейные зависимости, то среди чисел о появляются нули и соответствующие сингулярные векторы содержат коэффициенты этих зависимостей.
Используя разложение (10), можно решить прямоугольную (в общем случае несовместную) систему уравнений
J А« = АА(80).
(11)
Р, ] . (9)
Соотношение (9) дает обобщенное матричное представление всех частных случаев соотношений (6) и (7) в [37].
Программа в [37, 38], использующая процедуру масштабирования по Пулаи [24-31], основана в известной мере на программе Понгора [39]. В ней используется также соотношение (3) как критерий оптимума.
Рассмотрим более подробно метод минимизации, применяемый в [35-38]. Как известно, каноническая версия метода Гаусса-Ньютона, основанная на определении поправок А8 по соотношению (5), испытывает затруднения, если J почти сингулярна (почти вырождена): в этом случае поправки становятся неустойчивыми. Появление такой ситуации не означает вырождения задачи как таковой; она связана только со свойствами ее локальной линейной аппроксимации (4) в окрестности текущего приближения. Однако этого достаточно, чтобы понизить скорость сходимости или полностью нарушить сходимость итерационной процедуры. Эффективной процедурой для построения стабильного решения системы линейных
которая следует из (4). Приближенное решение для (11) отвечает минимуму суммы квадратов невязок. Это решение в то же время соответствует точному решению системы (5), поскольку (11) и (5) - не что иное, как условная и нормальная системы уравнений метода наименьших квадратов [41]. Заменяя I в (11) ее представлением (10), приходят к выражению вектора поправок А8 в виде суммы смещений х1 вдоль направлений, определяемых сингулярными векторами:
А« = V*. (12)
Величины х , соответствующие ненулевым сингулярным числам, связаны с элементами вектора ъ = иА^ соотношением х 1 = / ог Если о = 0, то соответствующую величину х1 можно выбрать произвольной. Следовательно, получают бесконечное множество решений системы с сингулярной матрицей. В частности, выбор х 1 = 0 дает решение А« с минимальной эвклидовой нормой (так как ортогональное преобразование в (12) сохраняет норму вектора). Эту ситуацию можно интерпретировать как исключение смещений вдоль сингулярных векторов, соответствующих нулевым о, при отыскании минимума. Такие смещения не изменяют величины линеаризованного функционала и, следовательно, они бесполезны.
Если матрица Якоби почти сингулярна, то система (5) формально имеет единственное решение. Однако
1
1
Р
Г
Г
эта матрица становится плохо обусловленной. В случае плохо обусловленной матрицы J все о отличны от нуля, но некоторые из них крайне малы по сравнению с другими. Замена малых сингулярных чисел нулями равносильна введению точных линейных соотношений между неизвестными вместо приближенных. Так как при этом некоторые базисные направления исключены, как и ранее, отыскание решения проводится в подпространстве более низкой размерности. Однако эта новая задача лучше обусловлена, чем первоначальная. Действительно, отношение наибольшего сингулярного числа к наименьшему ненулевому дает число обусловленности рассматриваемой задачи. Удаление минимального сингулярного числа понижает это отношение и, следовательно, приводит к улучшению обусловленности.
Применение в [35-38] указанного выше подхода обеспечивает устойчивую сходимость итерационного процесса Гаусса-Ньютона при нахождении масштабирующих множителей. Принцип отбора наименьших возможных поправок в случае их неоднозначного определения предотвращает появление расходимости.
При этом подчеркивается, что разложение (10) проводится заново на каждом шаге итерационной процедуры. Поэтому случайное вырождение матрицы Якоби в некоторых изолированных точках на ниспадающей траектории от начального приближения к минимуму будет сказываться на величинах и направлениях некоторых шагов, но это не изменит конечного решения, если оно хорошо локализовано. Только когда минимизируемый функционал имеет характер "оврага", а матрица J почти сингулярна на всех шагах в ходе итерационной процедуры (т.е. решение на самом деле не единственно или плохо обусловлено), видоизмененный метод Гаусса-Ньютона будет способствовать получению приемлемого результата, наименее отличающегося от заданного исходного приближения. (Это происходит потому, что выбраны наименьшие смещения из всех возможных на каждом итерационном шаге.)
Интересно отметить, что иная версия применения разложения по сингулярным числам в итерационной процедуре Гаусса-Ньютона предложена в [42]. Этот подход основан на введении демпфирующих множителей для соответствующих величин х а не на замене малых сингулярных чисел нулями. Следовательно,
когда полное вырождение отсутствует, сингулярные векторы, соответствующие приближенным линейным соотношениям, не исключены как возможные направления при поиске минимума, а лишь ограничены смещения вдоль этих направлений.
Авторы ссылок [10, 11] используют функционал (6) для определения масштабирующих множителей по методу Пулаи [24-31]. Конечно, штрафная функция, ответственная за близость искомого решения для силовых постоянных к исходному набору их квантовомеханических значений, при этом сохраняется в минимизируемом функционале. Для сопоставления результатов расчета частот колебаний, выполненных через посредство масштабирующих множителей, определенных функционалом (6) и с помощью функционала (3) [35-38], выбраны силовые поля гексаф-торэтана (С2Р6).
Численный пример определения масштабирующих множителей
Силовые поля молекулы С2Р6 рассчитаны в харт-ри-фоковском приближении (НР/6-3Ш*//НР/6-3Ш*) [1] и в приближениях, учитывающих электронную корреляцию. Последние используют второй порядок теории возмущений Мёллера-Плессе (МР2/6-3Ш*// МР2/6-3Ш*) [43] и теорию функционала плотности (В3ЬУР/6-3Ш*// Б31УР/6-3Ш*) [44-48].
Результаты, полученные в [10] с использованием функционала (6) [8-11], представлены в табл. 2. Эта таблица содержит также экспериментальные частоты колебаний молекулы С2Р6, взятые из [10, 49], и результаты расчетов, проведенных в [35, 36] при помощи функционала (3) и программы, упомянутой в [35, 36]. Расхождения между результатами расчета в приближении НР/6-3 Ш*//НР/6-3 Ш* в [10] и экспериментальными данными для двух частот вырожденных валентных колебаний формальных связей С-Р составляют 150 см-1 (у7 и у10 в табл. 2). Чрезмерно высокая невязка наблюдается в [10] также для у8. Эти примеры в особенности поучительны, так как метод масштабирования по Пулаи предназначен прежде всего для эмпирической коррекции квантово-механических силовых полей с тем, чтобы учесть эффекты электронной корреляции. Однако использование в [10] квантовомеханического силового поля, полученного с применением теории возмущений Мёллера-Плессе (принятие во внимание корреляционных эффектов) приводит к понижению различия
Таблиц а
Экспериментальные и рассчитанные частоты колебаний (см ') этаноподобной молекулы С2Ж6 с силовыми полями, полученными в приближениях НР/6-31С*//1Ш6-31С*, М Р2/6-31С * // М Р2/6-31С * и ВЗ ЬУ Р/6-31С*// ВЗ ЬУ Р/6-31С * и масштабированными по методу Пулаи с помощью программ, представленных
в [35, 36] и в [11]
V Отнесение Сим. Эксп.[10,49] НР/6-ЗЮ*//НР/6- ЗЮ* М Р2/6-3 1 О* //М Р2/6-ЗЮ* ВЗЬУР/б-ЗЮ*
[10] [35,36] [10] [35,36] [10]
1 "У(С-Р)вал. а и 1417 1425 1439 1423 1443 1420
2 V (С-С )вал. 807 774 786 778 787 794
3 8(СР3)деф. 348 336 338 341 344 346
4 Т(С-С)тор. а 1 „ 68 69 68 66 65 68
5 V (С-Р)вал. а 2 и 1116 1101 1111 1121 1127 1125
6 8(СР3)деф. 714 698 700 698 698 705
7 "У(С-Р)вал. ее 1250 1109 1267 1221 1264 1249
8 8(СР3)деф. 619 553 608 611 615 620
9 Р (С Рз)маят. 372 350 376 384 378 387
10 "У(С-Р)вал. е„ 1250.5 1093 1268 1207 1268 1232
11 8(СРз)деф. 522 489 516 526 514 529
12 Р (С Рз)маят. 219 186 206 213 213 216
Примечание. Результаты расчётов в приближении НР/6-ЗШ*//НР/6-ЗШ*: длины связей г(С-С) = 1.5264 А, г(С-¥) = 1.3110 А, углы С-С-Р = 109.79°°, Р-С-Р = 109.15°, Р-С-С-Р = 180.00° или 60.00°; полная энергия -672.384482 хартри; масштабирующие множители «(С-С)вал. = 0.7511, «(С-Р)вал. = 0.7767, х(СР3)деф. = 0.8415, х(СР3)маят. = 0.8061 и «(С-С)торс. = 0.9600. Результаты расчетов в приближении МР2/6-310*//МР2/'б-310*: длины связей г(С-С) = 1.5300 А, г(С-Р) = 1.3400 А, углы С-С-Р = 109.73°, Р-С-Р = 109.22°, Р-С-С-Р = 180.00° или 60.00°; полная энергия -673.652380 хартри; масштабирующие множители «(С-С)вал. = 0.8694, «(С-Р)вал. = 0.9382, х(СР3)деф. = 0.9999, х(СР3)маят. = 0.9992 и «(С-С)торс. = 0.9999.
Результаты расчетов в приближении В31ЛТ/6-ЗШ*//В31ЛТ/6-ЗШ*: длины связей г(С-С) = 1.5444 А, г(С-Р) = 1.3382 А, углы С-С-Р = 109.84°, Р-С-Р = 109.10°, Р-С-С-Р = 180.00° или 60.00°; полная энергия -675.254865 хартри; масштабирующие множители «(С-С)вал. = 0.9726, «(С-Р)вал. = 0.9631, х(СР3)деф. = 1.0502, х(СР3)маят. = 1.0365 и «(С-С)торс. = 1.9555.
между исходным и искомым силовыми полями. Вклад штрафной функции в функционал (6) соответственно уменьшается, и находятся более приемлемые значения масштабирующих множителей и рассчитанных частот колебаний У7 и у10 молекулы С2Б6 (табл. 2). Тем не менее отклонение У10 от экспериментального значения все еще составляет более 40 см 1 (см. табл. 2, а также табл. 19, стр. 279 в [10]).
Заметим, что применение программы, рассмотренной в [35, 36], для расчета масштабирующих множителей для силовых полей С2Б6 в приближениях
нр/6-3т*//нр/6-3т* и мР2/6-3Ш*//мР2/6-3Ш*,
где возможную плохую обусловленность задачи можно устранить посредством ограничения пути поиска минимума для функционала (3), а не видоизменением критерия оптимума, как в функционале (6), приводит к невязкам частот У7 и у10, которые составляют около 18 см-1. Невязка, полученная в случае частоты У8 еще ниже (см. табл. 2). Следует отметить, что размерность этой задачи довольно низка, она хорошо обусловлена, и в этом случае ограничения на выбор решения в программе [35, 36], возможно, не являются необходимыми.
В случае силового поля, полученного в приближении В3ЬУР/6-3Ю*//В3ЬУР/6-3Ю*, учет корреляционных эффектов производится эмпирическим методом [44-48]. Это приводит к недооценке большинства рассчитанных частот колебаний по сравнению с их экспериментальными аналогами. Величина частоты торсионного колебания т(С-С)(оге оказалась существенно заниженной (46 см-1 по сравнению с экспериментальным значением 68 см-1 [10, 49]). Расчет масштабирующих множителей в [10] с использованием функционала (6), несмотря на более близкие значения теоретических силовых постоянных к оптимальным, приводит к отклонению частоты валентного колебания У10 приблизительно на 20 см 1 (см. табл. 2). Отметим также, что различие между рассчитанными частотами существенно совпадающих частот колебаний у7 (её) и у10 (еи) в [10] также составляет около 20 см 1 Это ясно указывает на искажение колебательных мод вследствие масштабирования по критерию (6). Различие в частотах колебаний У7 (её) и У10 (еи) связано, очевидно, с использованием системы зависимых внутренних координат в [10]. При этом переход от Г в декартовых коорди-
д-т
натах к матрице в зависимых внутренних координа-
тах (Г сопровождается наложением дополни-
тельных ограничений, требующих, чтобы:
1) норма всех недиагональных силовых постоянных квадратной матрицы силовых постоянных была минимальной;
2) ранг этой матрицы также был минимальным (см., например, [5, 13, 49-51]).
Очевидно авторы [10] в качестве минимального ранга этой матрицы имеют в виду 3К-6 (или 3К-5) в зависимых внутренних координатах.
Вновь первое дополнительное ограничение с очевидностью следует из идеологии решения обратной колебательной задачи (см. предположение 1 в разделе "Традиционный подход"). Второе условие вполне естественно, так как после получения матрицы силовых постоянных с минимальной нормой для всех недиагональных силовых постоянных в зависимой системе координат ее норма увеличена. Следует отметить, что ранее предложены два теоретически обоснованных дополнительных условия для перехода к зависимой системе координат, а именно так называемое каноническое силовое поле Кучеры [52] и регулярное силовое поле Пупышева и др. [53]. Определенные преимущества использования регулярного силового поля в зависимых координатах отмечены в [53].
Обсуждение
При решении обратной колебательной задачи с Гд-т [5, 8-10, 49, 54] внимание сосредоточено на том, как минимально следует изменить Гд-т в терминах эвклидовой нормы матрицы, чтобы получить наилучшее соответствие экспериментальным частотам колебаний. Присутствие штрафной функции а||Г-Гд-т||2 в функционале (6), по существу, сводится к добавлению некоторых исходных параметров, которые должны быть воспроизведены искомой матрицей Г. В то же время использование наименьшей нормы а||Г-Гд-т||2 в качестве критерия не накладывает очень жестких ограничений на каждую отдельную силовую постоянную и приводит к приравниванию значимости отклонений как для диагональных, так и для недиагональных элементов этих матриц. Однако точность определения этих элементов зависит от кван-товомеханического приближения и качества использованного базисного набора. Следует помнить, что значения недиагональных силовых постоянных иногда на один или даже два порядка ниже, чем значе-
ния диагональных, и приравнивание значимости их абсолютных отклонений приводит прежде всего к большим относительным изменениям в недиагональных элементах. Отметим, что эта особенность противоречит представлению о таком преобразовании Г —т, чтобы норма всех недиагональных силовых постоянных была бы минимальной [5, 13, 4951] (см. выше). Таким образом, эта итерационная процедура должна приводить к некоторым существенным изменениям в конечной матрице силовых постоянных Г (как к повышению, так и к понижению значений и изменению знаков ее элементов) по сравнению с Г д-т. Следствием этих изменений являются заметные искажения колебательных мод Ь и, следовательно, рассогласованность с физическими требованиями к отнесениям частот колебаний (см. обсуждение в [23, 30]). Действительно, для бензола [13] решение колебательной задачи с использованием Г д-т и функционала (6а) дает отклонения от экспериментальных частот порядка 200 см-1. Большие отклонения получены также для межъядерных расстояний в [13] от геометрических параметров, найденных из вращательных постоянных [17]. Поэтому применение функционала (6), основанного на начальном кванто-вомеханическом наборе силовых постоянных, для решения обратной колебательной задачи характеризует постановку задач в [5, 8-10, 49, 54] как ошибочную с точки зрения физической модели.
Скорее всего, уточнять геометрические параметры молекулы и определять эффекты ангармоничности нельзя на теоретическом уровне МР2/6-3Ш*. Кроме того, одновременное варьирование геометрии при решении колебательной задачи с исходным квантовоме-ханическим силовым полем означает выход из минимума на квантовомеханической поверхности потенциальной энергии. Очевидно, поэтому в [13] получено сокращение кольца бензола при учете эффектов ангармоничности. Положение не улучшается в присутствии штрафной функции Р||Я - КпШа1||2, так как структура довольно далека от искомой структу-
ры (см. выше).
Пренебрежение малыми силовыми постоянными, соответствующими дальним "взаимодействиям" на языке МВСП, как это сделано в [5, 49], также приводит к изменению знаков некоторых недиагональных силовых постоянных и их относительных значений в ходе итерационной процедуры [49]. Однако
при подгонке рассчитываемых частот к экспериментальным используются производные силовых постоянных по частотам колебаний. Эти производные могут быть довольно большими для малых силовых постоянных, соответствующих дальним "взаимодействиям" в Г .
д-т
Действительно, существенные изменения в значениях распределения потенциальной энергии (которые представляют собой функцию от колебательных мод) благодаря использованию функционала (6) для решения обратной колебательной задачи показаны для У3 и У6 (см. отнесения в табл. 2) молекулы С2Б6 в табл. 7 в [49]. Следовательно, в этом случае применение традиционной идеологии решения обратной колебательной задачи к квантовомеханическому силовому полю также оказывается ошибочным [49]. Это, несомненно, справедливо также для рекомендаций по наложению дальнейших ограничений на кванто-вомеханическое силовое поле [55], требующих в частности, чтобы норма всех недиагональных силовых постоянных была минимальна [5, 13, 49-51]. Указанные ограничения [55] применяются для матрицы силовых постоянных в зависимых координатах, которая подвергается преобразованиям. Использование зависимых координат также связано с идеологией решения обратной колебательной задачи. Дело в том, что неудачи при переносе силовых постоянных, определенных решением обратной колебательной задачи, на родственные молекулы были приписаны использованию системы независимых координат. В настоящее время это соображение реализовано как прямое масштабирование индивидуальных силовых постоянных в полном наборе обычных естественных зависимых координат [56]. Действительно, некоторые элементы матрицы силовых постоянных в этих координатах представляют собой определенные линейные комбинации силовых постоянных в зависимых координатах. Поэтому представлялось целесообразным переносить на родственные молекулы силовые постоянные, полученные решением обратной колебательной задачи в зависимых координатах. Однако квантово-механические расчеты показали, что силовые постоянные родственных молекул могут заметно отличаться друг от друга. Ряд родственных молекул, таких как бута-1,3-диен и некоторые олигоены [57-62], могут служить примером. Отметим, что относительные ошибки квантовомеханических силовых посто-
янных приблизительно одинаковы для одинаковых функциональных групп в рядах родственных молекул. Это является причиной, по которой их квантово-механические силовые поля родственных молекул можно исправить одинаковыми масштабирующими множителями. Кроме того, значения этих масштабирующих множителей близки друг к другу и могут быть применены для коррекции элементов матрицы силовых постоянных, состоящих из линейных комбинаций силовых постоянных в зависимых координатах с тем, чтобы получить результаты с приемлемой спектроскопической точностью. Особенно показательным примером переноса масштабирующих множителей на родственные молекулы являются производные 3,3-диметилциклопропена, содержащие изова-лентные гетероатомы, т.е. атомы одной и той же группы периодической системы Менделеева [63-67].
Хорошо известно [22, 68, 69], что использование так называемого гомогенного масштабирования, т.е. одного и того же масштабирующего множителя для всех силовых постоянных, не меняет матрицу колебательных мод Ь. Применение масштабирования по методу Пулаи приводит к набору масштабирующих множителей с чрезвычайно малым разбросом от их среднего значения. Это означает, что процедура масштабирования практически не меняет исходные ко -лебательные моды, рассчитанные с Р—т [22, 68, 69]. Это обстоятельство свидетельствует о выполнении основного физического критерия для матрицы масштабированных силовых постоянных, т.е. сохранения отнесений частот колебаний, полученных при решении прямой колебательной задачи с Р ч_т.
Применение функционала (6) в методе Пулаи для определения масштабирующих множителей, когда их число в отличие от решения обратной колебательной задачи понижается до <п, должно, вообще говоря, гарантировать постоянство знаков силовых постоянных при масштабировании Р —т. Однако такое определение масштабирующих множителей неявно накладывает дополнительные ограничения на искомую матрицу Р. Поэтому присутствие второго члена в функционале (6) может привести к масштабирующим множителям, отличным от оптимальных с точки зрения минимума для выражения (3). Действительно, если использование Р —т дает отклонения рассчитанных частот от их экспериментальных аналогов на ~20%, то включение штрафной функции
сс||Р-2 в функционал (6) даже для малых коэффициентов а, полученных обычным методом [8-11] на каждом итерационном шаге, приводит к несовместной задаче вследствие этого набора исходных параметров. В самом деле, искомая матрица Р (см. соотношение (7)) должна быть весьма близкой к Р—т, что дает довольно большие отклонения рассчитанных частот от их экспериментальных аналогов. В то же время искомая матрица Р должна точно воспроизводить частоты экспериментального колебательного спектра. Вследствие этого масштабирующие множители в [10], рассчитанные в приближении НР/6-3т*//НР/6-3Ш* и МР2/6-3Ш*//МР2/6-3Ш*, для молекулы С2Р6 дают частоты У7 и у10, значения кото -рых существенно ниже, чем экспериментальные частоты (табл. 2). Когда теоретическое приближение расчета силового поля повышается, то в случае С2Р6 отклонения чисто теоретических частот колебаний от их экспериментальных аналогов понижаются, т.е. в задаче на масштабирующие множители это приводит к уменьшению невязок частот колебаний. Этот эффект в частности очевиден для высоко симметричной молекулы С2Р6, так как размерности симметричных блоков матрицы Р не превышают третьего порядка (см. табл. 2). Для большей молекулы, обладающей более низкой симметрией, размерности симметричных блоков матрицы Р много выше и число недиагональных силовых постоянных резко возрастает. Число масштабирующих множителей также возрастает. Поэтому для большой молекулы этот эффект использования штрафной функции а||Р - Р 9_т|| 2не столь очевиден. Это проявляется, например, в случае молекулы 4-фторбензальдегида [50], где отклонения чисто теоретических частот (рассчитанных с Р ) от их экспериментальных аналогов, как и в случае С2Р6, составляют 200 см-1. Однако расчеты масштабирующих множителей и частот молекулы 4-фтор-бензальдегида в одном из случаев применения функционала (6) приводят к приемлемым невязкам частот колебаний (см. Расчет I в [50]). Тем не менее это не означает, что применение функционала (6) не сказывается на значениях найденных масштабирующих множителей.
При решении обратной колебательной задачи с использованием функционала (6), где ограничения, неявно накладываемые при определении масштабирующих множителей, отсутствуют и недиагональные
силовые постоянные могут менять свои знаки и свои соотношения, указанная проблема невязок, как правило, не возникает. Если функционал (6) использован для решения обратной колебательной задачи для молекулы С2Р6 с силовым полем Г д-т, рассчитанным в приближении НР/6-3 Ш*//НР/6-3 Ш* в качестве исходного приближения, тогда удается приближенно воспроизвести экспериментальные значения частот колебаний.
Однако это получено при допущении как изменений в соотношениях между недиагональными и диагональными силовыми постоянными, так и изменений в знаках недиагональных силовых постоянных по сравнению с Гд-т [49]. Другими словами, при этом элементы матрицы Г теряют свою иерархию по сравнению с Гд-т, а колебательные моды меняются. Полученные при этом результаты будут существенно зависеть от качества базисного набора и приближения квантовомеханического расчета.
Можно также отметить, что для молекулы среднего размера, подобной молекуле 4-фторбензальдегида, более целесообразно перенести масштабирующие множители, найденные для квантовомеханических силовых полей малых молекул с аналогичными структурными группами, на матрицу Г д-т этой молекулы и предсказать частоты колебаний с помощью полученной матрицы Г, а не рассчитывать масштабирующие множители, используя весь набор экспериментальных частот колебаний этой молекулы с заданными отнесениями. Для молекулы среднего размера определение масштабирующих множителей последним методом, т.е. с заранее принятым отнесением частот колебаний, равносильно решению обратной колебательной задачи, когда частоты колебаний заранее определены и отнесены.
Уместно также дать некоторые общие замечания относительно использования приближенной теории функционала плотности (БРТ) [44-48] в вычислительной колебательной спектроскопии. Несмотря на относительно успешное применение БРТ, коррекция соответствующих теоретических силовых полей представляет интерес до настоящего времени. Действительно, при использовании локальной БРТ для расчета частот колебаний бензола [70] оказалось, что рассчитанные частоты у14 (В2и) легкого и тяжелого бензолов отклоняются от хорошо известных экспериментальных частот на 70 см 1 (4,35%) и 88 см 1
(6,8%) соответственно. Невязка в 90 см 1 получена в [71] для анилина. В случае бута-1,3-диена рассчитанные с помощью силового поля В3ЬУР/6-3 Ш*// В3ЬУР/6-3Ш* частоты валентных колебаний формальных связей С=С завышены на 100 см 1 [72]. К сожалению, при использовании БРТ нельзя предсказать, какие рассчитанные частоты колебаний могут дать большие отклонения от экспериментальных значений. Нельзя также предсказать знаки таких больших отклонений от экспериментальных частот коле -баний. Более того, в случае больших молекул, где плотность экспериментальных полос может быть довольно высокой в определенных областях спектра (например, большое количество полос в интервале 100 см х), возможные существенные отклонения рассчитанных частот не позволяют выполнить правильное отнесение. Эта особенность БРТ резко понижает ее предсказательные возможности для колебательных спектров сложных молекул. Поэтому для силовых полей в приближении БРТ также требуются дополнительные эмпирические поправки.
Выводы
Основное отличие традиционного подхода от подхода, использующего масштабирование квантовоме-ханического силового поля, состоит в способе определения элементов матрицы Г. При решении обратной колебательной задачи требование близости искомого решения Г к исходному квантовомеханическому набору силовых постоянных Гд-т согласно минимальной эвклидовой норме а|| Г-Г д-т||2 приводит к изменениям знаков недиагональных силовых постоянных случайным образом и к нарушению иерархии их абсолютных значений по сравнению с матрицей Г . Когда встречаются такие изменения и нарушения, то они обусловливают изменения в модах каждого отдельного колебания. Это может проявляться как ошибочные отнесения в конечных результатах проведенного колебательного анализа исследуемой молекулы (см. обсуждение в [22, 23, 28, 32]).
Специальный метод устранения возможной плохой обусловленности, использованный в [8-11], и идеология, применяемая для чисто эмпирического решения обратной колебательной задачи, не противоречат друг другу (см. предположение 2 в разделе "Традиционный подход"). Однако этот метод и различные дополнительные ограничения, накладываемые на Гд-т [5, 13, 49-51, 55, 73], изменяют исход-
ное квантовомеханическое силовое поле, которое используется для решения обратной колебательной задачи, и переводит задачу в целом на чисто эмпирический уровень.
Следовательно, если применение функционала (6) к какой-либо задаче дает согласно терминологии авторов публикаций [5, 9, 10] "псевдорешение", тогда результаты расчетов в [13], полученных для структурных параметров бензола (табл. 1), можно рассматривать как "псевдоструктуры". Этот вывод подтверждается также тем результатом, что учет ангармоничности в [13] вызывает сокращение бензольного кольца.
Решение обратной колебательной задачи с Р [5, 8-10, 13] можно также рассматривать как определение максимального числа масштабирующих множителей для матрицы Р ц_т, потому что при этом используются различные масштабирующие множители для ее диагональных и недиагональных элементов (см. обсуждение в [30]).
Основой идеологии масштабирования по методу Пулаи [24-31] является предположение о приблизительно одинаковых относительных ошибках, привносимых в результаты квантовомеханических расчетов силовых постоянных для рядов родственных молекул. Это предполагает переносимость масштабирующих множителей в рядах родственных молекул. При этом основным физическим критерием является близость для исследуемой молекулы колебательных мод, полученных с масштабированным силовым полем, к колебательным модам, определенным с не масштабированным квантовомеханическим силовым полем. Этот критерий должен, по существу, выполняться, когда разброс значений масштабирующих множителей минимален.
Пример молекулы С2Р6 (см. также [35,36]) показывает, что требование близости искомого решения к исходному приближению, заимствованное из традиционной идеологии чисто эмпирического подхода к решению обратной колебательной задачи и реализованное в [10] путем включения штрафной функции в минимизируемый функционал, приводит к тому, что полученные масштабирующие множители не могут должным образом исправить матрицу Рц-щ. Поэтому масштабированные силовые поля, полученные в [10], неверно воспроизводят экспериментальные частоты колебаний (см. табл. 2). Это означа-
ет, что перенос такого набора масштабирующих множителей на силовое поле родственной молекулы большего размера препятствует правильному предсказанию частот ее колебательного спектра. Очевидно, силовые постоянные, найденные таким образом, не будут переносимы в ряду родственных молекул.
Для малых высокосимметричных молекул, таких как С2Р6, результаты расчетов в [10] объясняются следующим образом. При решении обратной колебательной задачи условие близости к Р согласно эвклидовой норме не накладывает ограничений, которые могли бы привести к неразрешимости задачи. В качестве нежелательного результата при этом следует ожидать изменения знаков недиагональных элементов и нарушения иерархии их абсолютных значений в искомой матрице Р [49]. В случае ограничений, неявно наложенных при отыскании масштабирующих множителей, и когда Рц т дает отклонения квантовомеханических частот от экспериментальных колебательных частот на 200 см-1, включение штрафной функции а||Р-Рц_т||2 делает задачу практически несовместной и, следовательно, не решаемой. Когда качество базисного набора улучшено и теоретическое приближение при расчете матрицы Р повышено, расчет дает квантовомеханические частоты колебаний, которые ближе к экспериментальным значениям, и вклад эвклидовой нормы а||Р-Р ||2 понижается. Поэтому искомая матрица Р лучше воспроизводит экспериментальные частоты колебаний ( табл. 2).
Таким образом, конструкции, подобные функционалу (6), пригодны для определения масштабирующих множителей только в случае, когда Рц_т достаточно близка к искомой матрице Р. Напротив, процедура масштабирования по методу Пулаи предназначена для эмпирических поправок матрицы Рц_т, рассчитанной вблизи хартри-фоковского предела, на эффекты электронной корреляции [27]. Конечно, учет эффектов электронной корреляции такими методами, как теория возмущений Мёллера-Плессе (МР) и БРТ, могут уменьшить различие между матрицами Р и Р. Однако применение процедуры масштабирования по методу Пулаи в этом приближении будет довольно формальным, так как неизвестно, пригодны или нет такие поправки к силовым постоянным, для которых эффекты электронной корреляции частично учтены.
Заметим, что дополнительное условие близости конечной матрицы Г к Г д-т не следует из рассмотрения физической подоплеки метода масштабирования [24-31]. Параметр регуляризации а в (6) и (6а) не обладает также каким-либо физическим значением. Поэтому алгоритм регуляризации, приведенный в [8-11], адекватен требованиям чисто эмпирического решения обратной колебательной задачи и не пригоден для поправок квантовомеха-нического силового поля. Действительно, чисто механический перенос идеологии традиционного эмпирического решения обратной колебательной задачи путем использования функционала (6) на определение масштабирующих множителей и коррекцию квантовомеханического силового поля этими масштабирующими множителями неизбежно приводит к недоразумениям.
Отметим также, что задача определения масштабирующих множителей сама по себе, как правило, хорошо определена в отличие от обратной колебательной задачи (см., например, [74]). Тем не менее специальный алгоритм для проверки обусловленнос-
ти задачи должен быть включен в расчетную программу, так как на практике гарантирует получение устойчивого решения. Ввиду всех комментариев, данных для каждого из рассмотренных выше подходов, метод разложения по сингулярным числам матрицы Якоби [40, 41] представляется как наиболее приемлемая процедура для повышения устойчивости решения при отыскании масштабирующих множителей для коррекции матрицы Г
В общем, традиционная идеология вычислительной колебательной спектроскопии основана на интуитивных моделях силовых полей и на решении обратной колебательной задачи, в то время как современная идеология вычислительной колебательной спектроскопии базируется на неэмпирическом квантовомеханическом расчетном методе, который развит из первых принципов, и на эмпирических поправках к соответствующим результатам. Если первая идеология пытается навязать свои жесткие модели природе, то последняя, опирающаяся на познанные законы природы, способствует дальнейшему ее изучению.
СПИСОК ЛИТЕРАТУРЫ
1. Hehre W.J., Radom L., Schleyer v. R. P., Pople J.A. II Ab Initio
Molecular Orbital Theory. Wiley, N.Y., 1986.
2. Волькенштейн М.В., Ельяшевич М.А., Степанов Б.И. II
Колебания молекул. M., 1949.
3. Вильсон Е., Дешиус Дж., Кросс П. II Теория колебательных
спектров молекул. M., 1960.
4. Сивин С. II Колебания молекул и среднеквадратичные
амплитуды. M., 1971.
5. Kuramshina G.M., Weinhold F., Pentin Yu.A. II J. Chem. Phys.
1998. 109. P. 7286.
6. Бард Й. II Нелинейное оценивание параметров. M., 1979.
С. 101.
7. Химмельблау Д. II Прикладное нелинейное программирова-
ние. M., 1975. С. 233.
8. Ha T.-K., Meyer R., GünthardHs.H. II Chem. Phys. Lett. 1978.
59. P. 17.
9. Кочиков И.В., Курамшина Г.М., Пентин Ю.А., Ягола А.Г. II
Обратные задачи колебательной спектроскопии. M., 1993.
10. Yagola A.G., Kochikov I.V., Kuramshina G.M., Pentin Yu.A. II Inverse Problems of Vibrational Spectroscopy. VSP BV. Utrecht. The Netherlands. 1999.
11. Кочиков И.В. , Курамшина Г. М., Степанова А.В., Ягола А.В. II Вестн. Моск. ун-та. Сер. 3. Физика. Астрономия. 1997. № 5. C. 21.
12. Дашевский В.Г., Баранов А.П., Медведь Т.Я., КабачникМ.И.
II ДАН СССР 1983. 270. № 2. С. 355.
13. Kochikov I.V., Tarasov Yu.I., Kuramshina G.M., Spiridonov V.P., Yagola A.G., Strand T.G. II J. Mol. Struct. 1998. 445. P. 243.
14. Gundersen S., Strand T.G, Volden H.V. II J. Mol. Struct. 1998. 445. P. 35.
15. Goodman L., Ozkabak A.G., Thakur S.N. II J.Phys.Chem. 1991. 95. P. 9044.
16. Huber K.P., Herzberg G. II Molecular Spectra and Molecular
Structure. IV. Constants of Diatomic Molecules. Van Nostrand Reinhold Company. N.Y., 1979.
17. Gauss J., Stanton J.F. II J. Phys. Chem. A. 2000. 104. P. 2865.
18. Махнев А.С., Степанов Н.Ф., Панченко Ю.Н. II Вестн. Моск. ун-та, Сер. 2. Химия. 1977. № 6. С. 644.
19. Махнев А.С., Степанов Н.Ф., Панченко Ю.Н., Нипан M.E., Maтвеев В.^Н Вестн. Моск. ун-та. Сер. 2. Химия. 1978. № 1. С. 16.
20. Ermakov K.V., Butayev B.S., Spiridonov V.P. II J. Mol. Struct.
1990. 240. P. 295.
21. Kuchitsu K., Nakata M., Yamamoto S. II Stereochemical Applications of Gas-Phase Electron Diffraction, in: I. Hargittai and M. Hargittai (Eds.) VCH. N.Y., Part A. 1988. P. 227.
22. Панченко Ю.Н., Степанов Н.Ф. II Журн. физ. химии. 1995. 69. C. 592.
23. Panchenko Yu.N. II J.Mol.Struct. 1997. 410-411. P. 327.
24. Panchenko Yu.N., Pulay P., TörökF. II J. Mol. Struct. 1976. 34. P. 283.
25. Fogarasi G., Pulay P. II Annu. Rev. Phys. Chem. 1984. 35. P.
191.
26. Fogarasi G., Pulay P. II Vibrational Spectra and Structure I Ed. J.R. Durig.Elsevier. Amsterdam. 1985. 14. P. 125.
27. Pupyshev V.I., Panchenko Yu.N., Bock Ch. W., Pongor G. II J. Chem. Phys. 1991. 94. P. 1247.
28. Panchenko Yu.N., De Mare G.R., Stepanov N.F. // J. Mol. Struct.
1995. 348. P. 413.
29. Panchenko Yu.N., De Mare G.R., Pupyshev V.I. // J. Phys. Chem.
1995. 99. P. 17544.
30. Панченко Ю.Н. // Изв. РАН. 1996. № 4. C. 800.
31. Panchenko Yu.N.// J. Mol. Struct. 2001. 567-568. P. 217.
32. Stepanov N.F., De Mare G.R., Panchenko Yu.N. // J. Mol. Struct. (Theochem). 1995. 342. P. 9.
33. Гельман Г.Г. // Квантовая химия. М.Л.; 1937. Гл. IX. C. 428.
34. Feynman R.P. // Phys. Rev. 1939. 56. P. 340.
35. Panchenko Yu.N., De Mare G.R., Stepanov N.F. // Russ. J. Phys.Chem. 2000. 74. Suppl. 2. P. 245.
36. Панченко Ю.Н., Абраменков А.В. // ЖФХ. 2003. 77. C. 1062.
37. Краснощеков С.В., Абраменков А.В., Панченко Ю.Н. // Вестн. Моск. ун-та. Сер. 2. Химия. 1985. 26. № 1. C. 29.
38. Краснощеков С.В., Абраменков А.В., Панченко Ю.Н. // ЖФХ.
1997. 71. C. 497.
39. Pongor G. // Program SCALE. Eötvös Lorand Univ., Dept. of General and Inorg. Chem., Budapest, Hungary, 1978, and its two new versions.
40. Уилкинсон Дж., Райнш К. // Справочник алгоритмов на языке Алгол. Линейная алгебра. M., 1976.
41. Лоусон Ч., Хенсон Р. // Численное решение задач метода наименьших квадратов. M., 1986.
42. Brodersen S. // J. Mol. Spectr. 1990. 142. P. 122.
43. Moller C, PlessetM.S. // Phys. Rev. 1934. 46. P. 618.
44. KryachkoE.S., LudenaE.V. // Energy Density Functional Theory of Many-Electron System. Kluwer. 1990.
45. Joubert D. (ed.) // Density Functionals: Theory and Applications.
Springer, 1998.
46. Ziegler T. // Chem. Rev. 1991. 91. P. 651.
47. Becke A.D., Jr. // J. Chem. Phys. 1993. 98. P. 5648.
48. Lee A.D., Yang W., Parr R.G. // Phys. Rev. 1988. B37. P. 785.
49. Курамшина Г.М., Вэйнхолд Ф., Кочиков И.В., Пентин Ю.А., Ягола А.В. // ЖФХ. 1994. 68. № 3. C. 401.
50. Kochikov I.V., Tarasov Yu.I., Spiridonov V.P., Kuramshina G.M., Yagola A.G., Saakjan A.S., Popik M.V., Samdal S. // J. Mol. Struct. 1999. 485-486. P. 421.
51. Kochikov I.V., Tarasov Yu.I., Spiridonov V.P., Kuramshina G.M., Saakjan A.S., Yagola A.G. // J. Mol. Struct. 2000. 550-551.
P. 429.
52. Kuczera K. // J. Mol. Struct. 1987. 160. P. 159, ibid, idem, 1984.
117. P. 11.
53. Pupyshev V.I., Panchenko Yu.N., De Mare G.R., Bock Ch. W. // J. Mol. Struct. 1992. 272. P. 145.
54. Пентин Ю. А., Курамшина Г. М. // Журн. структ. хим. 1995. 36. C. 204.
55. Курамшина Г. М., Ягола А.Г. // Журн. структ. хим. 1997. 38. C. 221.
56. Baker J., Jarzecki A.A., Pulay P.// J. Phys. Chem. A. 1998. 102. P. 1412.
57. Bock Ch. W., Panchenko Yu.N., Krasnoshchiokov S. V., Pupyshev
V.I. // J. Mol. Struct. 1985. 129. P. 57.
58. Panchenko Yu.N., Krasnoshchiokov S. V., George Ph., Bock Ch. W. // Struct. Chem. 1992. 3. P. 15.
59. Panchenko Yu.N., Bock Ch.W. // Struct. Chem. 1992. 3. P. 27.
60. Panchenko Yu.N., De Mare G.R., Aroca R., Bock Ch. W. // Struct Chem. 2000. 11. P. 121.
61. Panchenko Yu.N., Vander Auwera J., Moussaoui Ya., De Mare G.R. // Struct Chem. 2003. 14. P. 337.
62. De Mare G.R., Panchenko Yu.N., Vander Auwera J. // J. Phys. Chem. A. 1997. 101. P. 3998.
63. Panchenko Yu.N., De Mare G.R. // Spectrochim. Acta. 2003. 59A. P. 329.
64. Panchenko Yu.N., De Mare G.R., Abramenkov A. V., Baird M.S.,
Tverezovsky V.V., Nizovtsev A.V., Bolesov I.G. // Spectrochim. Acta A. 2003. 59А. P. 2087.
65. De Mare G.R., Panchenko Yu.N., Abramenkov A. V., Baird M.S., Tverezovsky V.V., Nizovtsev A.V., Bolesov I.G. // Spectrochim. Acta A. 2003. 59А. P. 2063.
66. Panchenko Yu.N., De Mare G.R., Abramenkov A. V., Baird M.S.,
Tverezovsky V.V., Nizovtsev A.V., Bolesov I.G. // Spectrochim. Acta A. 2003. 59А. P. 1733.
67. De Mare G.R., Panchenko Yu.N., Abramenkov A. V., Baird M.S.,
Tverezovsky V.V., Nizovtsev A.V., Bolesov I.G. // Spectrochim. Acta A. 2004. 60А. P. 519.
68. Pupyshev V.I., Stepanov N.F., Krasnoshchiokov S.V.,
De Mare G.R., Panchenko Yu.N. // J. Mol. Struct. 1996. 376. P. 363.
69. Краснощеков С.В., Степанов Н.Ф., Панченко Ю.Н. // Журн. структ. хим. 1998. 39. C. 210.
70. Berces A., Ziegler T. // J. Chem. Rev. 1993. 98. P. 4793.
71. Rauhut G., Pulay P. // J. Phys. Chem. 1995. 99. P. 3093.
72. Zhou X., Wheelees C.J.M., Liu R.. // Vibr. Spectrosc. 1996. 12. P. 53.
73. Kuramshina G.M., WeinholdF. // J. Mol. Struct. 1997. 410-411. P. 457.
74. Кочиков И.В., Курамшина Г.М., Ягола А.Г. // Журн. вычисл. и матем. физ. 1987. 27. C. 1651.
Поступила в редакцию 28.06. 04
SCALING THE QUANTUM MECHANICAL MOLECULAR FORCE FIELD VERSUS SOLVING THE INVERSE VIBRATIONAL PROBLEM Yu.N. Panchenko, G.R. De Mare
(Division of Physical Chemistry)
The peculiarities characterising the traditional approach (solving the inverse vibrational problem) used in calculational vibrational spectroscopy and the unconventional approach based on using scaled quantum mechanical molecular force fields are considered. Some results on the determination of the equilibrium geometry of benzene in both the harmonic approximation and in the approximation taking into account the kinematic and dynamic anharmonicity corrections by solving the inverse vibrational problem are discussed. Using the quantum mechanical force fields of the C2F6 molecule, calculated at three different theoretical levels as an example, the results of the determination of scale factors by different mathematical techniques are compared.
