χ' = S-1/2χ ; (3.9.3)
или, что то же самое, вместо матрицы DS используется матрица S1/2DS1/2.
Методы Малликена и Лёвдина дают разные заряды на атомах, но математически оба они равноценны – нельзя сказать, какой из них даёт лучшие результаты. У метода Малликена больше всего недостатков, основные из которых следующие:
1) диагональные элементы матрицы DS могут быть больше 2, что противоречит принципу Паули;
2) недиагональные элементы могут быть отрицательными, что невозможно с физической точки зрения;
3) деление недиагональных вкладов поровну между атомами ничем не оправдано; более логично, например, делить эти вклады в соответствии с электроотрицательностью атомов;
4) базисная функция может иметь малую экспоненту, что соответствует хорошему описанию электронной плотности далеко от этого атома, а электронная плотность на этой функции полностью приписывается данному атому.
Значения атомных зарядов зависят от метода их вычисления, поэтому с их помощью нельзя рассчитать дипольные, квадрупольные и другие моменты.
У метода Лёвдина первые три из вышеперечисленных недостатков отсутствуют.
Матрица плотности может быть также использована для расчёта порядков связей. В соответствии с моделью Малликена сила связи между атомами А и В может быть определена из матрицы плотности:
Порядок связи по Малликену (; 3.9.4) также зависит от базисного набора, но в меньшей степени, чем заряды на атомах.
3.10. Анализ заселённостей на основе электростатического потенциала
Значительная часть несвязных взаимодействий между полярными молекулами описывается как электростатическое взаимодействие между фрагментами с несимметрично распределённой электронной плотностью. С фундаментальной точки зрения рассматривается взаимодействие между электростатическим потенциалом (ESP), генерированным одной молекулой, и заряженными частицами другой. Электростатический потенциал в данной точке вычисляется как сумма вкладов ядер и электронной волновой функции:
Первая часть потенциала вычисляется обычным способом из зарядов ядер и их координат. Электронный вклад требует знания волновой функции.
Электростатический потенциал напрямую зависит от волновой функции. Поэтому при вычислениях наблюдается его сходимость по мере увеличения размера базиса и степени вычисления энергии корреляции. Так как потенциал напрямую зависит от электронной плотности, то в общем он не сильно зависит от уровня вычислений. Например, на уровне Хартри – Фока в расщеплённом поляризованном базисе DZP получаются достаточно хорошие результаты.
Стоуном (Stone) был разработан метод распределённого мультипольного анализа (DMA – Distributed Multipole Analysis). Этот метод основан на том факте, что электростатический потенциал, возникающий из зарядового перекрывания двух базисных функций, может быть описан в виде ряда мультиполей относительно некоторой точки между двумя ядрами. Эти моменты можно рассчитать непосредственно из матрицы плотности и базисных функций. Мультипольный ряд конечен. Наибольший неисчезающий член определяется суммой орбитальных моментов двух функций. Например, две р-функции дают мультиполи только до квадруполя. Для гауссовых орбиталей центр мультиполей между ядрами А и В вычисляется из координат ядер и гауссовых экспонент по простой формуле:
Если подобные вычисление мультиполей провести для каждой пары функций, то мы получим точный электростатический потенциал. Таких вычислений нужно провести М2. Для элементов 1-го и 2-го периодов достаточно рассмотреть s- и р-функции. При этом центры зарядов, диполей и квадруполей будут располагаться на ядрах или на центрах связей.
3.11. Анализ заселённостей на основе
собственно волновой функции
Волновая функция позволяет рассчитать вероятность нахождения электрона в определённой точке пространства, т.е. волновая функция даёт пространственное распределение заряда. Для вычисления зарядов на атомах в молекуле остаётся только разделить всё её пространство на части, принадлежащие отдельным "атомам". Но атомов фактически в молекулах нет, т.к. квантовая механика рассматривает молекулу как совокупность ядер и электронов. Поэтому нужно решить, как разделить пространство молекулы на части, соответствующие разным ядрам. При этом нет необходимости рассматривать базисные функции, как в методе Малликена или Лёвдина.
Метод разделения пространства молекулы на атомные области был предложен Бэйдером (Bader) и получил название метода "атомов в молекуле" (AIM – Atoms In Molecule). В этом методе атомные области разделены граничными поверхностями, для которых градиент электронной плотности по отношению к ближайшим ядрам равен нулю, т.е. на этой поверхности электрон испытывает одинаковое притяжение к обоим ближайшим ядрам. Таким образом формируются области вокруг ядер, которые считаются атомами в составе молекулы. Далее проводится интегрирование квадрата волновой функции по соответствующим областям и находятся соответствующие заряды.
Ещё один метод анализа волновой функции предложен Циословским (Cioslowski). Полученные этим методом заряды на атомах называются зарядами обобщённого атомного полярного тензора (GAPT – Generalized Atomic Polar Tensor).
Атомный полярный тензор (APT – Atomic Polar Tensor) представляет собой матрицу 3 3 производных дипольного момента по координатам (3.11.1). Он определяет интенсивность поглощения инфракрасного излучения.
Эта матрица зависит от выбора системы координат, но её след не зависит от этого. Поэтому Циословски предложил вычислять заряд на атоме А как 1/3 следа этой матрицы (3.11.2):
Так как дипольный момент сам по себе является первой производной от энергии, то необходимо вычислять соответствующие вторые производные. Это требует большого объема вычислений. Но если проводится расчёт ИК-спектров, для чего необходимо вычислять эти вторые производные (матрицу Гессе), то попутно вычисляются и соответствующие заряды.
Кроме больших вычислительных затрат, недостатком зарядов GAPT является их чувствительность к степени учёта корреляции. Эти два недостатка ограничивают использование зарядов GAPT.
3.12. Локализованные орбитали
Волновую функцию Хартри – Фока можно написать в виде детерминанта Слэйтера, составленного из набора ортонормированных МО. Для расчётных целей лучше использовать канонические МО, т.е. те, которые приводят матрицу множителей Лагранжа к диагональному виду и которые являются собственными функциями оператора Фока. Но после завершения процедуры самосогласования мы можем выбрать другой набор орбиталей путём составления линейных комбинаций канонических МО (3.12.1). Общая волновая функция и все наблюдаемые свойства будут идентичны для такого преобразования (вращения) МО.
Традиционно считается, что химическая связь объясняется возрастанием вероятности нахождения электрона между двумя ядрами в сравнении с суммой вкладов чистых атомных орбиталей. Канонические МО делокализованы по всей молекуле и не отражают этого. Канонические МО не отражают также и концепцию функциональных групп. Одной из целей вычисления локализованных молекулярных орбиталей (ЛМО) является получение таких МО, которые будут почти неизменны для подобных структурных (функциональных) групп в различных молекулах. Набор ЛМО может быть получен оптимизацией значения математического ожидания двухэлектронного оператора (3.12.2):
Практически локализацию проводят с помощью ряда орбитальных вращений размерности 2 × 2. Так как наблюдаемые свойства зависят от общей электронной плотности, а не от плотности на некоторой МО, то двухэлектронный оператор можно выбрать различными способами.
В схеме локализации Бойса (Boys) в качестве двухэлектронного оператора используется квадрат расстояния между двумя электронами, т.е. набор ЛМО определяется путём минимизации пространственных экспонент (3.12.3):
Схема Эдмистона – Руденберга (Edmiston – Ruedenberg) использует в качестве двухэлектронного оператора величину, обратную расстоянию между двумя электронами (3.12.4), что соответствует определению набора ЛМО путём нахождения максимума энергии отталкивания:
Метод фон-Ниссена (von Niessen) использует дельта-функцию от расстояния между двумя электронами в качестве оператора и максимизирует его математическое ожидание:
Метод Пайпека – Мезей (Pipek-Mezey) соответствует нахождению максимума суммы атомных зарядов Малликена, т.е. находят максимум следующей функции:
Реже всего используется схема фон-Ниссена. Остальные три метода дают очень близкие результаты. Главным исключением являются системы, содержащие как σ-, так и π-связи, – например этилен. Схема Пайпека – Мезей сохраняет различия между σ- и π-связями, в то время как схемы Бойса и Эдмистона – Руденберга дают изогнутые "банановые" связи. С другой стороны, для плоских молекул с неподеленными парами, например для воды, схема Пайпека – Мезей даёт одну неподеленную пару в плоскости молекулы, а вторую – вне плоскости π-типа. Методы Бойса и Эдмистона – Руденберга дают в этом случае правильную картину – обе неподеленных пары эквивалентны и представляют собой "заячьи уши".
Скорость вычислений методом Бойса пропорциональна , для метода Эдмистона – Руденберга – , для метода фон-Ниссена – также , для метода Пайпека – Мезей – .
В методе Эдмистона – Руденберга используется более фундаментальный энергетический критерий локализации по сравнению с другими методами. В методе Бойса критерием локализации служат расстояния, в методе Пайпека – Мезей – заряды. Но методы Бойса и Пайпека – Мезей более быстрые и поэтому они используются чаще, тем более что полученные разными методами ЛМО мало отличаются.
Натуральные орбитали получают диагонализацией матрицы плотности. Собственные векторы матрицы плотности представляют собой натуральные орбитали, а собственные значения – орбитальные заселённости. Для функции Хартри – Фока заселённости могут быть равны 0 или 2, для корреляционных функций – любому значению от 0 до 2. Анализ заселённостей на основе натуральных орбиталей проводится аналогично процедуре Лёвдина. При этом общая волновая функция строится в виде детерминанта (или детерминантов) Слэйтера, составленных из натуральных атомных орбиталей.
3.13. Примеры: метан и вода
Из табл. 3.1 и 3.2 видно, что методы Малликена и Лёвдина не приводят к какому-либо определённому значению для зарядов с увеличением базиса. Значения зарядов в этом случае непредсказуемы. Если в базис включены диффузные функции, то результат вообще может быть абсурдным.
Таблица 3.1
Атомные заряды на углероде в СН4 [4]
Базис
|
Малликен
|
Лёвдин
|
ESP
|
NAO
|
AIM
|
STO-3G
3-21G
6-31G(d,p)
6-311G(2d,2p)
6-311++G(2d,2p)
cc-pVDZ
cc-pVTZ
cc-pVQZ
aug-cc-pVDZ
aug-cc-pVTZ
|
-0.26
-0.80
-0.47
-0.14
-0.18
-0.13
-0.37
-0.27
-0.63
-1.20
|
-0.15
-0.38
-0.43
-0.13
-0.20
-0.76
-0.21
-0.07
-0.43
+0.05
|
-0.38
-0.45
-0.36
-0.36
-0.36
-0.31
-0.35
-0.36
-0.35
-0.37
|
-0.21
-0.89
-0.88
-0.69
-0.71
-0.79
-0.72
–
-0.77
-0.72
|
+0.25
-0.01
+0.26
+0.19
+0.19
+0.32
–
–
+0.33
–
|
Таблица 3.2
Атомные заряды на кислороде в Н2О [4]
Базис
|
Малликен
|
Лёвдин
|
ESP
|
NAO
|
AIM
|
STO-3G
3-21G
6-31G(d,p)
6-311G(2d,2p)
6-311++G(2d,2p)
cc-pVDZ
cc-pVTZ
cc-pVQZ
aug-cc-pVDZ
aug-cc-pVTZ
|
-0.39
-0.74
-0.67
-0.52
-0.47
-0.29
-0.48
-0.51
-0.26
-0.41
|
-0.27
-0.46
-0.44
-0.00
-0.12
-0.58
-0.11
+0.23
-0.39
+0.12
|
-0.65
-0.90
-0.81
-0.74
-0.76
-0.76
-0.75
-0.75
-0.74
-0.74
|
-0.41
-0.87
-0.97
-0.91
-0.93
-0.91
-0.92
–
-0.96
-0.93
|
-0.89
-0.93
-1.24
-1.24
-1.25
-1.27
–
–
-1.26
–
|
Заряды, вычисленные методами электростатического потенциала (ESP), "атомов в молекуле" (AIM) и в базисе натуральных АО, являются вполне определёнными величинами при увеличении базисного набора. Но, к сожалению, заряды, определённые этими тремя методами, значительно различаются между собой. Например, атому углерода в метане может быть приписан заряд от +0.2 (AIM) до -0.7 (NAO).
Таким образом, понятие заряда на атоме является теоретически неопределенным, и обсуждать заряды если и можно, то только для близких по структуре молекул, расчет которых проводился одинаковыми методами.
Таблица 3.3
Порядки связей в СН4 и Н2О [4]
Базис
|
CH4, DS
|
CH4, AIM
|
H2O, DS
|
H2O, AIM
|
STO-3G
3-21G
6-31G(d,p)
6-311G(2d,2p)
6-311++G(2d,2p)
cc-pVDZ
cc-pVTZ
cc-pVQZ
aug-cc-pVDZ
aug-cc-pVTZ
|
0.99
0.93
0.98
1.00
1.00
1.00
0.98
0.98
0.91
0.86
|
1.00
1.00
1.00
1.00
1.00
0.99
–
–
0.99
–
|
0.95
0.83
0.89
0.99
0.96
1.03
1.01
1.01
1.11
1.04
|
0.81
0.80
0.62
0.63
0.62
0.60
–
–
0.62
–
|
3.14. Полуэмпирические методы
Полуэмпирические методы предполагают последующее введение ряда дополнительных приближений в метод Хартри – Фока. Первым шагом является рассмотрение только валентных электронов. Роль óстовных электронов сводится к понижению эффективного заряда ядра. Более того, выбирается минимальный базис, позволяющий разместить валентные электроны. В результате водороды имеют по одной базисной функции, а атомы 2-го и 3-го периодов имеют только по 4 базисных функции – одну s- и три p-функции. Большинство современных полуэмпирических методов использует только s- и p-функции, которые выбираются в виде орбиталей Слэйтера, т.е. являются экспоненциальными.
Центральным звеном всех полуэмпирических методов является приближение нулевого дифференциального перекрывания (НДП, ZDO – Zero Differential Overlap). В этом приближении отбрасываются все произведения базисных функций, зависящие от координат одного и того же электрона и локализованные на разных атомах. Если обозначить атомную орбиталь на атоме А как μа, то это приближение математически можно выразить как μа(i) . νb(i) = 0. Необходимо отметить, что нулю считается равным произведение функций, а не интеграл этого произведения. Это приближение приводит к следующим упрощениям:
1) матрица перекрывания S становится единичной;
2) одноэлектронные трёхцентровые интегралы (два центра от базисных функций и один от оператора) становятся равными нулю;
3) все трёх- и четырёхцентровые двухэлектронные интегралы также считаются равными нулю и отбрасываются.
Чтобы компенсировать эти три приближения, остальные интегралы считаются параметрами и их значения вычисляются на основе экспериментальных данных или приравниваются к каким-либо экспериментальным величинам. Таким образом, различия между полуэмпирическими методами определяются следующим: 1) сколько интегралов считается равными нулю; 2) каким образом проводится параметризация оставшихся интегралов.
В настоящее время наиболее распространен модифицированный метод пренебрежения двухатомным перекрыванием (Modified Neglect of Diatomic Overlap – MNDO), разработанный группой Дьюара (M.J.S. Dewar). В модифицированных методах частично отказываются от приближения нулевого дифференциального перекрывания, но интегралы вычисляют с помощью какой-либо параметризации. Модификации метода Дьюара различаются: 1) способом расчёта межъядерного отталкивания; 2) тем, как определяются параметры.
Одноэлектронные одноцентровые интегралы имеют значение, соответствующее энергии притяжения электрона к ядру (Us, Up) плюс вклады притяжения к остальным ядрам системы. Последние параметризуются в виде произведения зарядов (пониженных) ядер на двухэлектронный интеграл:
где и – s- или p-функции атома А.
Двухцентровые одноэлектронные интегралы вычисляются как произведение соответствующего интеграла перекрывания на среднее значение двух атомных "резонансных параметров" β:
где и – s- или p-функции атомов А и В; интеграл перекрывания рассчитывается численным интегрированием. Последнее не относится к приближению НДП, и поэтому метод называется "модифицированным".
В методе MNDO остаётся только 5 типов одноцентровых двухэлектронных интегралов в sp-базисе:
Параметры G-типа являются кулоновскими интегралами, а параметры Н – обменными интегралами; интеграл Gp2 включает две разных р-функции.
Межъядерное (остов-остовное) отталкивание в методе MNDO вычисляется по формуле (3.14.9):
где экспоненты α являются параметрами.
Взаимодействия, включающие O-H и N-H связи, рассчитываются по другой формуле (3.14.10):
В дополнение к этому MNDO предполагает равенство экспонент для s- и р-орбиталей: ζs = ζр. Параметры Gss , Gsp , Gpp , Gp2 , Hsp взяты из атомных спектров, а значения других параметров подобраны так, чтобы наилучшим образом воспроизводить свойства молекул.
Параметры метода MNDO оптимизированы для следующих элементов: H, B, C, N, O, F, Al, Si, P, S, Cl, Zn, Ge, Br, Sn, I, Hg и Pb.
Остиновская модель 1 (AM1 – Austin Model 1). После проведения расчётов методом MNDO стало ясно, что он даёт ряд систематических ошибок. Например, отталкивание между двумя атомами, находящимися на расстоянии 2-3 Å, слишком велико. Из этого следует, что энергии активации также слишком высоки. Это является следствием завышенного отталкивания остова. Поэтому óстовная функция была модифицирована добавлением гауссовых функций и вся модель была репараметризована. В результате получилась остиновская модель AM1, названная в честь переезда Дьюара в университет Остина. Остов-остовное отталкивание в AM1 имеет следующий вид (3.14.11):
где величина k равна от 2 до 4, в зависимости от типа атома.
Как и для MNDO, параметры Gss , Gsp , Gpp , Gp2 , Hsp взяты из атомных спектров. Другие параметры, включая константы ak , bk и ck, подобраны для наилучшего воспроизведения свойств молекул.
Параметры метода АM1 оптимизированы для следующих элементов: H, B, C, N, O, F, Al, Si, P, S, Cl, Zn, Ge, Br, Sn, I и Hg.
Параметризация методов MNDO и АМ1 была проведена вручную, вследствие чего было взято очень небольшое число соединений в качестве “обучающего” набора. Стюарт разработал способ автоматической параметризации, основанный на вычислении производных от ошибок по отношению к параметрам. Далее все параметры были вновь автоматически подобраны, включая двухэлектронные вклады. Обучающий набор был выбран достаточно большим – несколько сотен соединений. Для этой репараметризации формула метода АМ1 для остов-остовного отталкивания была сохранена, за исключением того, что для каждого атома оставлялись только две гауссовых функции (k = 2). В результате метод получил соответствующее название, которое далее было сокращено до РМ3 (MNDO-PM3 – MNDO, Parametric Method №3). По сути это тот же метод АМ1 с полной оптимизацией всех параметров. Можно считать, что получился наилучший набор параметров. Но на результат оптимизации повлиял и человеческий фактор: были отобраны для обучающего набора только определённые соединения, выбирался фактор значимости для каждого набора экспериментальных свойств.
Параметры метода РM3 оптимизированы для следующих элементов: H, Li, C, N, O, F, Mg, Al, Si, P, S, Cl, Zn, Ga, Ge, As, Se, Br, Cd, In, Sn, Sb, Te, I, Hg, Tl, Pb, Bi, Po, At.
Можно сформулировать следующие общие ограничения для методов MNDO, AM1, PM3 и им подобных.
1. Барьеры вращения вокруг связей, имеющих частично двоесвязанный характер, существенно занижены. Это касается вращения вокруг связи C-N в амидах, где для соответствующих барьеров получены значения всего 5-10 ккал/моль. Неэмпирический расчёт даёт величины 20-25 ккал/моль, что соответствует экспериментальным данным. Для барьера вращения вокруг центральной связи бутадиена полуэмпирические методы также дают низкое значение – только 0.5-2.0 ккал/моль, тогда как экспериментальное значение – 5.9 ккал/моль.
2. Слабые взаимодействия, такие как ван-дер-ваальсовы, водородные связи, предсказываются плохо. При этом и взаимодействие слишком слабое, и геометрия неправильная.
3. Слишком низкое значение длины связи с нитрозильной группой. Например, связь N-N в N2O3 получается короче на 0.7 Å.
4. Хотя методы MNDO, AM1 и PM3 имеют параметры для некоторых металлов, эти параметры вычислены лишь на небольшом наборе экспериментальных данных, поэтому нужно осторожно относиться к расчётам структур, содержащих атомы металлов.
Так как метод АМ1 содержит больше варьируемых параметров, чем MNDO, а метод РМ3 является версией АМ1 с полностью оптимизированными всеми параметрами, то следует ожидать, что ошибка будет уменьшаться в ряду:
MNDO > АМ1 > РМ3.
Но эта закономерность относится, во-первых, лишь к средней ошибке, и в каждом конкретном случае любой из методов может быть лучшим. Во-вторых, эта закономерность относится лишь к тем свойствам, по которым проводилась параметризация. Например, при вычислении зарядов на атомах метод АМ1 даёт более реальные результаты, особенно для соединений с атомами азота.
Для химиков наиболее интересны относительные энергии для ряда структур. Поэтому, несмотря на сравнительно хорошую точность воспроизведения теплот образования – 5-10 ккал/моль – это, к сожалению, лишь средняя ошибка. При сравнении энергий двух соединений средние ошибки должны складываться, что даёт уже более существенную ошибку – до 20 ккал/моль. Это является существенным отличием полуэмпирических методов от неэмпирических. Методы ab initio обычно лучше воспроизводят относительные энергии, чем абсолютные. Это объясняется преобладанием систематических ошибок в этих методах, которые компенсируются при сравнении энергий подобных систем.
4. ТЕОРИЯ ФУНКЦИОНАЛА ПЛОТНОСТИ
(Density Functional Theory)
Методы теории функционала электронной плотности (ДФТ) разработаны сравнительно недавно и также позволяют учесть электронную корреляцию. Но эти методы можно отнести к корреляционным лишь частично, так как они принципиально отличаются от всех методов квантовой химии.
Теория функционала плотности основана на доказательстве Хоэнбергом и Коном (Hohenberg, Kohn) следующей теоремы: электронная энергия основного состояния полностью определяется распределением электронной плотности ρ(r). Другими словами, существует однозначное соответствие между энергией и электронной плотностью:
E(r) ↔ ρ(r). (4.1)
Важность этого заключения можно показать из сравнения с анализом системы на основе волновой функции. Волновая функция для N-электронной системы содержит 3N координат – по 3 на каждый электрон. Электронная плотность является квадратом волновой функции, интегрированной по N-1 электронным координатам. Поэтому электронная плотность является функцией только трех координат независимо от числа электронов. Единственной проблемой остаётся то, что функционал, связывающий энергию и электронную плотность, неизвестен. Задача методов ДФТ состоит в построении подобного функционала.
Возможность применения метода ДФТ в компьютерной химии появилась с введением орбиталей Коном и Шамом (Kohn and Sham). Ключевым положением теории Кона – Шама является расчёт кинетической энергии в предположении невзаимодействующих электронов. Это является аналогией одноэлектронному приближению Хартри – Фока. В действительности электроны взаимодействуют между собой, и уравнение для функционала кинетической энергии TS не включает всю кинетическую энергию. Однако, аналогично методу Хартри – Фока, она учитывает 99% энергии. Разница между вычисленной и реальной энергией очень мала. Оставшуюся кинетическую энергию Кон и Шам включили в обменно-корреляционный член EXC[ρ]. В результате получилось следующее выражение для полной энергии в методе ДФТ:
Приравнивая EDFT к точной энергии, это выражение можно рассматривать как определение для обменно-корреляционной энергии. Это часть энергии, которая остаётся после вычитания кинетической энергии, не учитывающей взаимодействия электронов, а также энергии притяжения электронов к ядрам (Ene) и энергии кулоновского отталкивания (J):
Достоинством метода ДФТ является то, что необходимо рассчитать только полную электронную плотность. Но чтобы рассчитать кинетическую энергию, всё же необходимо ввести орбитали. Несмотря на это метод ДФТ по затратам времени сравним с методом Хартри – Фока, но может давать более точные результаты.
Главной проблемой в методе ДФТ является выбор формул для обменно-корреляционного вклада. Если функционал найден, то вычисления проводятся аналогично волновой механике: определяется набор ортогональных орбиталей путём минимизации энергии. Фактически проводится расчёт псевдо-собственных значений и собственных функций, называемых каноническими орбиталями Кона – Шама. Соответствующие уравнения называются уравнениями Кона – Шама:
где
Множители Лагранжа () здесь также можно рассматривать как энергии молекулярных орбиталей, а энергию высшей занятой орбитали считать энергией ионизации в соответствии с теоремой Купманса. Это справедливо только в случае точного обменно-корреляционного функционала, но так как в реальных вычислениях это невозможно, то орбитальные энергии имеют не совсем то же значение, как и в методе Хартри – Фока. Неизвестные орбитали Кона – Шама могут быть определены численными методами или распространены на набор базисных функций аналогично методу Хатри – Фока. Различия в методах ДФТ состоят в выборе формы функционала обменно-корреляционной энергии.
4.1. Типы функционалов – LDA, GGA и гибридные
В функционале приближения локальной плотности (LDA – Local Density Approximation) считается, что локально электронную плотность можно рассматривать как однородный электронный газ. Это эквивалентно тому, что плотность является медленно изменяющейся функцией. Обобщением метода LDA является метод LSDA, который отдельно рассматривает электронную плотность для электронов с разным значением спиновой переменной (LSDA – Local Spin Density Approximation).
Приближение LSDA в общем переоценивает обменную энергию примерно на 10%, что превышает ошибку, связанную с пренебрежением корреляции электронов. Электронная корреляция также переоценивается примерно в 2 раза, сила связи также завышается. Методы LSDA таким образом дают результаты, сравнимые по точности с методом Хартри – Фока.
Методы градиентной корректировки (GGA). Для усовершенствования метода LSDA необходимо рассматривать модель неоднородного электронного газа. Шагом в этом направлении является введение зависимости обменной и корреляционной энергии не только от электронной плотности, но и от производных от плотности. Подобные методы известны как методы градиентной корректировки или приближения обобщённого градиента (GGA – Generalized Gradient Approximation). Большинство из них связано с модификацией функционала LSDA. Методы получили названия по именам разработчиков и по году опубликования. Наиболее известные из них – PW86 и PW91 (Perdew, Wang), B88 (Becke), LYP (Lee, Yang, Parr). Последний функционал обычно используется для расчёта органических соединений. В методах градиентной корректировки в функционал плотности корреляционной энергии включены 4 параметра. Их значения подобраны так, чтобы наилучшим образом воспроизводились экспериментальные данные для атома гелия.
Наиболее совершенными в настоящее время являются гибридные методы. Из уравнения Гамильтона и определения обменно-корреляционной энергии можно получить точное соотношение между обменно-корреляционной энергией и соответствующим потенциалом, связывающим невзаимодействующую систему с реальной системой. Полученное таким образом уравнение называется формулой адиабатической связи, в которой проводится интегрирование по параметру λ, включающему электрон-электронное взаимодействие:
В первом приближении можно считать, что потенциал Vxc линейно зависит от λ и может быть представлен как среднее от двух предельных значений:
При λ = 0 получаем предельный случай, соответствующий невзаимодействующим электронам. При этом учитывается только обменная энергия, энергия корреляции отбрасывается. Точной волновой функцией в этом случае является один детерминант Слэйтера, составленный из орбиталей Кона – Шама. Обменная энергия точно соответствует теории Хартри – Фока. Величина второго вклада (соответствующего λ = 1) неизвестна. Применение данного подхода в рамках приближения LSDA называется методом "половина-на-половину" – Half-and-Half ("H+H"):
Модели, включающие точную обменную энергию, называются гибридными методами, которые соответствуют адиабатически связанной модели. Примером является функционал Бека, включающий 3 параметра:
Параметры a, b и c подбирают таким образом, чтобы лучше воспроизводить экспериментальные данные. Значения параметров определяются формой корреляционного функционала , и типичные их значения: . Процедура В3 была обобщена путём включения большего числа параметров, что, однако, не дало существенных улучшений.
4.2. Преимущества и недостатки методов ДФТ
По точности результатов методы GGA сравнимы с методом МР2 или лучше него. Но методы GGA требуют меньше времени счёта – столько же, сколько обычный метод Хартри – Фока. Метод ДФТ позволяет провести вычисления для систем, к которым метод МР2 неприменим, а полученные результаты при этом сравнимы с результатами кластерных методов.
Другим достоинством является то, что методы ДФТ, основанные на неограниченных детерминантах, не подвержены спиновой засорённости. Это является следствием того, что электронная корреляция оказывается включённой в однодетерминантную функцию (в виде Ехс). Другим следствием включения электронной корреляции в однодетерминантную функцию является устойчивость вычислительного процесса к потере симметрии неограниченного детерминанта по сравнению с функциями Хартри – Фока. Например, при расчёте молекулы озона нельзя получить более низкую энергию, соответствующую функции UHF для "чистых" методов ДФТ (LSDA, BLYP, BPW91). Но функционалы, включающие точную обменную энергию (B3LYP, B3PW91), показывают триплетную нестабильность молекулы озона.
Слабые взаимодействия (например, ван-дер-ваальсовы) плохо описываются существующими функционалами. Например, функционал LDA переоценивает силу связи и предсказывает притяжение между атомами разреженного газа. Водородные связи имеют в основном электростатическую природу и хорошо описываются методами ДФТ. По некоторым данным, относительные энергии не очень хорошо предсказываются методами ДФТ и переходные состояния плохо описываются с помощью этих методов.
И, наконец, методы ДФТ плохо описывают возбуждённые состояния, имеющие одинаковую с основным состоянием симметрию. Отсутствие волновой функции в описании системы делает затруднительным гарантировать ортогональность между основным и возбуждённым состояниями.
Величины абсолютных ошибок методов ДФТ минимальны для функционала B3LYP, который используется наиболее часто (табл. 4.1).
В общем можно считать, что методы ДФТ, особенно градиентные и гибридные, являются значительно более точными по сравнению с полуэмпирическими группы MNDO. В методах ДФТ ошибки ближе к систематическим. И, таким образом, методы ДФТ можно считать хорошим инструментом для расчёта систем, где высокая точность не требуется.
Таблица 4.1
Достарыңызбен бөлісу: |