Вычислительные методы

Для систем без сферической симметрии явное построение функции отклика независимых частиц (2.24) численно требовательно, так как она является функцией двух координат гиг'. Это особенно актуально в методах сеток, таких как разложения в базисе плоской волны и представления координатной сетки. По этой причине были разработаны некоторые эффективные вычислительные методы, которые позволяют избежать явного построения функции отклика [38]

2.2.1 Представление 3D Декартовой сетки

В методе конечных разностей в 3D Декартовых координатах координаты задаются дискретными равномерно. Обычно рассматриваются только валентные электроны, а электрон-ионное взаимодействие заменяется псевдопотенциалом, сохраняющим норму. Уравнения Кона-Шэма в дискретном виде могут быть получены с дискретными Декартовыми координатами х_} = lh, y_m = mh, z_n = nh, где h - шаг сетки, al, m, n целые числа:

h

• - Л^С1',т;,п'(Ш +l'h,y_m+m’h,z_n + n'h)
• 2m i

+^vioc(^xi>y_m>^znW(^xi>y_m>^zn) . (2.26)

+ Z^Vnonloc (^X1 ’ Ут ’ ^zn ! ^ХГ > Ут'> ^zn' Ш^X1'’ Ут'’ ^zn')

Г,пГ,п’

^{= ?}$(^х1>Ут>^гп)

Локальный потенциал V_loc включает потенциал Хартри, обменно-корреляционный потенциал и локальную часть ионного потенциала. Нелокальный потенциал V_nonIoc представлен нелокальной частью псевдопотенциала.

Чтобы уменьшить число точек сетки при сохранении высокой численной точности, следует использовать формулы конечных разностей высокого порядка для оператора кинетической энергии, обычным значением является 9-точечная формула для каждой Декартовой координаты.

2.2.2 Эволюция волновой функции в реальном времени

Гамильтониан Кона-Шэма зависит от времени, но на небольшом промежутке времени эту временную зависимость можно игнорировать. Поэтому эволюция во времени от t до t-ь At может быть вычислены в момент времени t из Г амильтониана:

• (t + At) = exp(-ih_KS(t)At /Й)^,(t) .
• (2.27)

После разложения в ряд Тейлора:

if

^i(t) •

(2.28)

Стабильность алгоритма зависит от k_max и шага времени At. Например, алгоритм й² Г л-у является стабильным при k_m;tY =4, при этом At < 2.87//^_nviY. где ~ 3 v ITlclA " г П IclA " IlldX 1

2т у h у

максимальное собственное значение Гамильтониана h_KS [38].

2.2.3 Метод реального времени

Можно рассчитать поляризуемость, решив уравнение нестационарные уравнения Кона-Шэма в режиме реального времени. Рассмотрим индуцированную поляризуемость р для электрического поля в произвольный момент в направлении v, V_ext (r,t) = eE(t)r_v. В случае линейного отклика каждый компонент w в p_A(t) и E(t) должен быть связан с дипольной поляризуемостью a_/7V(w):

jdtexp(iwt)p_A(t) = a_//H(w)jdtexp(iwt)E(t) . (2.29)

Поэтому, для данного приложенного электрического поля можно вычислить дипольную поляризуемость произведя преобразования Фурье:

«_/;r(^w)= g^fdtexp(iwt)p_/;(t) . (2.30)

где E(w) = Jdtexp(iwt)E(t) - преобразование Фурье приложенного электрического ПОЛЯ.

Есть два простых варианта для временного профиля электрического поля.

Первым способом является импульсное электрическое поле. Потенциал импульсного электрического поля:

V_ext(r,t) = l?(t)r_v , (2.31)

где I - величина импульса. В классической механике все электроны получают скорость от этого импульса, в квантовой механике каждая орбиталь умножается на плоскую волну с импульсом I. Перед получением импульса система находится в основном состоянии, которое описывается статическими орбиталями Кона-Шэма. Орбитали после получения импульса связаны с орбиталями до получения импульса соотношением:

• 0+

i//,(t₀₊) = exp

fdt(h_KS[no] +I^(t)r_v)^i(t₀_) .

(2.32)

Л о-

Для небольшого промежутка времени Гамильтониан Кона-Шэма может быть проигнорирован. Таким образом, орбитали после получения импульса:

(r,t = 0₊) = exp(ikr_v)^ (г) , (2.33)

где k = I / Й.

Начиная с этой начальной волновой функции орбитали развиваются без какого-либо внешнего потенциала.

Дипольный момент рассчитывается из орбиталей:

^d_A4) = fdrr^XhCrd)!² • (2.34)

i

Поляризуемость вычисляется как:

2 ⁰⁰

a_Ap(w) = - jdtexp(iwt)d_A(t) .

• (2.35)
• 1 о

Так как импульсное поле включает в себя равномерно все компоненты частоты, поэтому можно получить поляризуемость для всей спектральной области через эволюцию по времени.

Другим способом является применение однородного статического поля при t > 0:

E(t) =

0(t < 0)

E_o(t>O) ’

(2.36)

Преобразование Фурье задается как E(w) = -E₀/iw. До t = 0 система находится в основном состоянии, которое описывается статическими орбиталями Кона-Шэма. После включения электрического поля поляризация осциллирует около статической поляризации р^, которая связана со статической поляризуемостью

^{= б2}//г(^)Е() • Поляризуемость может быть рассчитана как:

^0) = “ Jdtexp(iwt)p_/Z(t) .

(2.37)

^Ьо о

Переход к связанным возбужденным состояниям представляется как дельта функция в мнимой части поляризуемости. На практике, время эволюции производится в течение ограниченного периода Т . Разрешение энергии спектра ЛЕ, полученного из преобразования Фурье, связано с периодом как ЛЕ ~ h/T .

В поляризации p_/z(t) переходы связанных возбужденных состояний представляется в виде осцилляций, которые сохраняются без затухания. Преобразование Фурье подобных колебаний в ограниченном периоде сопровождается покачиванием вокруг энергии возбуждения из-за резкого прекращения колебаний. Для того чтобы исключить подобное, добавляют дополнительную функцию затухания в преобразование Фурье. Это можно сделать добавив небольшую мнимую часть i/ в частоту w, что равносильно умножению на функцию затухания в преобразование Фурье в (2.30):

f(t) = exp(-i_Z) . (2.38)

Можно использовать другую функцию затухания, например:

Г t у Г t у

f(t) = l-3- + 2- . (2.39)

Эта функция обладает следующими свойствами: f(0) = l,

f’(0)= f(T)= f’(T) = O. Такое условие f(T) = f’(T) = O делает поляризацию p_A(t) плавно исчезающей при t=T. Условие f’(0) = 0 гарантирует, что интеграл силы осциллятора не изменится после введения функции затухания. Если использовать данную функцию затухания, то полная ширина на уровне половинной амплитуды (FWHM) находится как:

^fwhm ~6.2Й/Т . (2.40)

Рисунок 2.1 - Распределение силы осциллятора фуллерена С60 [17]. Использовалось адиабатическое приближение локальной плотности (ALDA)

С учетом (2.39) можно записать уравнение для поляризуемости (2.35):

е^2Т a_AV(w) = -—jdtf(t)exp(iwt)d_A(t) .

(2.41)

Сила осциллятора:

2w

S(w)= —Ima(w) .

(2.42)

Эта функция соответствует спектру оптического поглощения.

На рисунке 2.1 изображена зависимость силы осциллятора рассчитанная методом реального времени для фуллерена С60 [17].

На рисунках 2.2-2.5 изображены результаты расчетов методом реального времени для фуллеренов С60, С20, С80 в программном комплексе Octopus 4.01. Параметры расчета: сетка имеет вид сферы радиусом 1.5 нм, шаг времени 0.0007 hl эВ, расстояние между точками в сетке 0.003 нм. На рисунке 2.2 изображена зависимость дипольного момента (2.34) после воздействия импульса электрического поля для фуллерена С60. Изображено распределение силы осциллятора (2.42) для фуллерена С60 (рис. 2.3), С20 (рис. 2.4), С80 (рис. 2.5).

• 6.00 ->
• 4.00 -
• 2.00 - I

.—. 0 ( У 1.0 2.0 3.0 4.0 5.0 6.0 7* 8.0 9Т) * 0.

•ii- -2.00 ? •

.0.0

• -2.00
• -4.00 ?

t, hbar/eV

• -6.00
• -8.00
• -10.00

Рисунок 2.2 - Зависимость изменения дипольного момента фуллерена С60 от времени после воздействия импульса электрического поля

800.00

Рисунок 2.3 - Распределение силы осциллятора фуллерена С60, рассчитанное методом реального времени в программном комплексе Octopus

Рисунок 2.4 - Распределение силы осциллятора фуллерена С20, рассчитанное методом реального времени в программном комплексе Octopus

Есть несколько характерных особенностей метода реального времени. Во-первых, не нужно рассчитывать незанятые орбитали Кона-Шэма. Информация о незанятых орбиталях приходит через эволюцию по времени орбиталей Кона-Шэма. Так как рассчитывается эволюция во времени для занятых орбиталей, то требования к памяти в расчетах равно умноженному на два количеству памяти для расчета статических орбиталей Кона-Шэма. Во-вторых, можно получить спектр всей частотной области от единого времени эволюции. В-третьих, вычислительный алгоритм достаточно прост [17, 23, 37, 38].

1400.0 и

Рисунок 2.5 - Распределение силы осциллятора фуллерена С80, рассчитанное методом реального времени в программном комплексе Octopus

2.2.4 Метод реального времени для периодических систем

Оптические свойства объемных периодических систем характеризуется частотно-зависимой диэлектрической функцией.

Для периодических систем выражения, полученные для конечных систем, не могут применены, поэтому сначала нужно сделать калибровочное преобразование:

ie tic

(2.43)

где A(t) - нестационарный пространственно-однородный векторный потенциал, связанный с электрическим полем как:

E(t) = -(l/c)dA(t)/dt .

(2.44)

Нестационарное уравнение Кона-Шэма (TDKS) с внешним электрическим полем E(t) преобразуется в:

А² ~A(t)f_p С )

^{+ V}loe №i(0

(2.45)

+ exp(-ieA(t)r_r /ftc)V_nonloc exp(ieA(t)r_r /йс)^ (t)

где f_v - единичный вектор в направлении v.

Так как дипольный момент не может быть определен для бесконечных систем, поэтому вместо него рассматривается ток интегрированный по пространству, которое является производной по времени дипольного момента:

• 1Д‘) =
• 7~РД0 = - ^е

dt m ; ^J dtr„

• (2.46)
• 1 ^Yr1/
• -i |Zf drdr>* (г,1)У_поп1ос(г’Г ^r//Ynonloc (r,r')^(r,t)

где нелокальная часть псевдопотенциала вносит свой вклад в ток, как указано в последнем члене.

Для применения метода реального времени для бесконечных периодических систем необходимо принимать во внимание эффект поверхностного заряда. Макроскопическое электрическое поле внутри материала E(t) состоит из внешнего электрического поля, которое пропорционально электрической индукции D(t), и макроскопического электрического поля, вызванного поляризацией P(t) на поверхности:

E(t) = D(t)-4/zP(t) . (2.47)

Так как поляризация на поверхности индуцируется током внутри элементарной ячейки, то:

lp(t) = li(t), (2.48)

dt Q ^Л

где О. представляет собой объем элементарной ячейки, a i_/z(t) определяется как интеграл тока по элементарной ячейке объемом О.. Пусть индуцированный векторный потенциал Д_П(1 (t), который связан с поляризацией:

p(t)= ¹ dAnd(t) ₍₂₄₉₎

4лс dt

Тогда уравнение (2.48) преобразуется:

д² . 4яс. ,_ч

(2.50)

~ 2 ^^nd ^) ^— гл ’ dr Q

Векторный потенциал A(t) в нестационарном уравнении Кона-Шэма (2.45) представляет собой сумму индуцированного векторного потенциала Д_пс1 (t) и внешнего векторного потенциала 4_xt(t), который связан с электрической индукцией D(t),_Ka_K D(t) = -(l/c)(54_xl(t)/5t).

Как и при расчете поляризуемости, можно вычислить диэлектрическую проницаемость ?(w) из расчета уравнений Кона-Шэма в реальном времени с внешним векторным потенциалом. Диэлектрическая проницаемость ?(w) определяется из соотношения:

jdtexp(iwt)D(t) = ^(w)jdtexp(iwt)E(t) .

(2.51)

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

• 1
• ?(w)

j dt exp(i wt)

(2.52)

J dt exp(iwt) где A_ot(t) = 4_xt(t) + 4_ld(t).

Внешний векторный потенциал можно задать через ступенчатую функцию 4xt (О = 4^(0 > ^что соответствует импульсному электрическому полю E(t) = —(l/c)(dA(t)/dt) = -(l/c)4