Примеры

5.3.1 Исследование физических свойств молекулы водорода

Состоит из нескольких частей (файлы tbasel- tbase5), которые позволяют рассчитать полную энергию, длину связи, плотность заряда и энергию атомизации.

На рисунке 5.2 представлен код, который позволяет рассчитать полную энергии молекулы водорода. После выполнения создается выходной файл tbasel_l.out.

Вычисление оптимального расстояния между атомами молекулы может быть выполнено несколькими способами:

• 1) вычислить энергию на разных расстояниях между атомами и найти такое расстояние, где энергия минимальна;
• 2) вычислить силы на разных расстояниях между атомами и найти такое расстояние, где сила равна нулю;
• 3) произвести процедуру минимизации энергии (или найти точку, когда сила равна нулю).
• ?Определение различных наборов данных ndtset 21

xcart: -0.5 0.0 0.0 # Стартовое значение атомных координат

0.5 0.00.0

xcart+ -0.025 0.0 0.0 # Приращение к xcart

0.025 0.0 0.0 # getwfk -1 nband 1

acell 10 10 10

ntypat 1 znucl 1

natom 2 typat 1 1

ecut 10.0

kptopt 0 nkpt 1

nstep 10 toldfe 1.0d-6 diemac 2.0 optforces 1

Рисунок 5.3 - Конфигурационный файл для расчета энергии на разных расстояниях между атомами водорода (tbase 1_2.in)

На рисунке 5.3 представлен код для расчета энергии на разных расстояниях между атомами водорода. При выполнении данного скрипта происходит вычисление суммарной энергии в зависимости от межатомного расстояния с шагом 0.05 радиуса Бора. В выходном файле в результате вычислений получится список с величинами энергий на разных межатомных расстояний, минимальное значение будет соответствовать оптимальному расстоянию.

На рисунке 5.4 представлен код. в котором производится структурная оптимизация положений атомов в молекуле водорода. При выполнении данного скрипта происходит автоматическое нахождения минимума с пределом сходимости силы для остановки вычисления. В выходном файле после нахождения минимума будут указаны координаты атомов в этой точке.

ionmov 3 # Структурная оптимизация

nt iine 10

tolmxf 5.0d-4 xcart -0.7 0.0 0.0

0.7 0.0 0.0

toldff 5.0d-5 nband 1 acell 10 10 10

ntypat 1 znucl 1

natom 2 typat 1 1

ecut 10.0

kptopt 0 nk.pt 1

nstep 10

diemac 2.0

Рисунок 5.4 - Конфигурационный файл для расчета структурной оптимизации молекулы водорода (tbasel_3.in)

На рисунке 5.5 представлен код для расчета энергии изолированного атома водорода. Ранее была вычислена суммарная энергия молекулы. Здесь вычисляется энергия одного изолированного атома водорода. Энергия атомизации будет находиться как разница между энергией молекулы и удвоенной атома водорода [1].

nsppol 2

occopt 2

nband 1 1

occ 1.0 0.0

spinat 0.0 0.0 1.0

acell 10 10 10

ntypat 1 znucl 1

natom 1

typat 1 xcart

0.0 0.0 0.0

ecut 10.0

kptopt 0

nkpt 1

nstep 10

toldfe 1.0d-6 diemac 2.0

optforces 1

Рисунок 5.5 - Конфигурационный файл для расчета энергии изолированного атома водорода (tbasel_5.in)

5.3.2 Нахождение конвергентных значений длины связи и энергии атомизации для молекулы водорода

В предыдущем примере были получены значения для длины связи молекулы водорода (1.522 Боровских радиуса) и энергии атомизации (4.506 эВ), которые сильно отличаются от экспериментальных (1.401 Боровских радиуса и 4.747 эВ).

На рисунке 5.6 представлен код, при выполнении которого производится расчет оптимальной длины связи и энергии атомизации в зависимости от параметра ecut (отсечка кинетической энергии). Были рассмотрены значения параметров ecut от 10 до 35 Хартри.

Другой способ - изменение размеров ячейки (tbase2_3.in), при этом производятся расчеты оптимальной длины связи и энергии атомизации в зависимости от параметра acell, сторона ячейки изменялась от 8 до 18 Боровских радиуса. Оптимальными были значения ecut=30 Хартри и сторона ячейки равная 12 107

Боровским радиусам. Дальнейшее увеличение этих параметров не вело к увеличению точности.

ndtset 12 udtset 6 2 acell 10 10 10

ecut:? 10 ecut+? 5

?Первый набор данных: нахождение оптимальной дли связи Н2 и полной энергии

?Второй набор данных: получение полной энергии отдельного атома natom?2 1 nsppol?2 2 occopt?2 2 nband?2 1 1 осс?2 1.0 0.0

toldfe?2 1.0d-6 xcart?2 0.0 0.0 0.0 spinat?2 0.0 0.0 1.0

ntypat 1 znucl 1 typat 1 1 kptopt 0 nkpt 1 nstep 10 diemac 2.0 optforces 1

Рисунок 5.6 - Конфигурационный файл для расчета оптимальной длины связи и энергии атомизации молекулы водорода (tbase2_2.in)

Для увеличения точности расчетов можно использовать приближения обобщенного градиента (GGA) вместо приближения локальной плотности (LDA). Для этого нужно поменять псевдопотенциалы и параметр ixc (например, ixc=ll).

В результате были получены значения для длины связи, которые отличались на 2-4 % от экспериментальных, а для энергии атомизации - на 2-5 % [1].

5.3.3 Исследование физических свойств кристаллического кремния

• ?определение типов атомов ntypat 1 znucl 14
• ?определение атомов natom. 2

typat 1 1 xred

• 0.0 0.0 0.0
• 1/4 1/4 1/4
• ?определение набора данных плоских волн ecut 8.0
• ?определение сетки к-точек

kptopt 1

ngkpt 222

nshiftk 4

?определение процедуры сходимости nstep 10

toldfe 1.0d-6

diemac 12.0

optforces 1

Рисунок 5.7 - Конфигурационный файл для нахождения полной энергии кристаллического кремния (tbase3_l.in)

В данном примере рассчитывается полная энергия, параметры решетки и зонная структура кристаллического кремния. На рисунке 5.7 представлено содержание конфигурационного файла для нахождения полной энергии кристаллического кремния. Для улучшения сходимости нужно использовать сетку с другим количеством точек (оптимальным значением ngkpt является «6 6 6»).

ndtset 2

?Набор данных 1 : самосогласованный расчет kptoptl 1 nshiftkl 4 shiftkl 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5

ngkptl 444

prtdenl 1 toldfel 1.0d-6

?Набор данных 2 : расчет зонной структуры iscf2 —2 getden2 -1 kptopt2 -3 nband2 8

ntypat 1 znucl 14 natom 2 typat 1 1 xred

• 0.0 0.0 0.0
• 1/4 1/4 1/4

ecut 8.0 nstep 10

diemac 12.0 optforces 1

Рисунок 5.8 - Конфигурационный файл для расчета зонной структуры кристаллического кремния (tbase3_5.in) [1]

На рисунке 5.8 представлено содержание конфигурационного файла для расчета зонной структуры кристаллического кремния. Зонная структура вычисляется путем решения уравнений Кона-Шэма для различных k-точек по различным направлениям зоны Бриллюэна. В результате получено значение для ширины запрещенной зоны 0.594 эВ, что далеко от экспериментального значения 1.17 эВ. Для того, чтобы получить правильные значения нужно использовать приближение GW.

5.3.4 Исследование физических свойств металла

occopt 4 tsmear 0.05

?задание ячейки acell 3*7.60

rprim 0.0 0.5 0.5

• 0.5 0.0 0.5
• 0.5 0.5 0.0
• ?оптимизация параметров решетки optcell 1 ionmov 3 ntime 10 dilatmx 1.05 ecutsm 0.5
• ?определение типов атомов ntypat 1 znucl 13
• ?определение атомов

natom 1

typat 1

xred

• 0.0 0.0 0.0
• ?определение набора данных плоских волн ecut 6.0
• ?Обменно-корреляционный функционал ixc 1
• ?определение сетки к-точек ngkpt 222 nshiftk 4

?определение процедуры сходимости nstep 10

toldfe 1.0d-6

Рисунок 5.9 - Конфигурационный файл для оптимизация параметров решетки кристаллического алюминия (tbase4_l.in)

В данном примере рассчитывается полная энергия, параметры решетки, для кристаллического алюминия.

На рисунке 5.9 представлено содержание конфигурационного файла, при выполнении которого производится оптимизация параметров решетки кристаллического алюминия.

Далее производится исследование сходимости. Одним из вариантов улучшить сходимость является увеличение числа к - точек (параметр ngkpt до значений «8 8 8»). В другом случае предлагается изменять число к - точек и одновременно параметр tsmear (значения от 0.01 до 0.04 Хартри) [1].

5.3.5 Расчет корректировки собственной энергии для собственных значений Кона-Шэма в приближении GW

В примере 5.33 была рассчитана зонная структура кремния, но при этом для ширины запрещенной зоны был получен далекий от эксперимента результат. Для того, чтобы более точно вычислить ширину запрещенной зоны, используется приближение GW.

Рисунок 5.10 - Блок-схема расчета в приближении GW

Расчет в приближении GW выполняется за четыре шага (рисунок 5.10):

• 1. Нахождение электронной плотности в основном состоянии.
• 2. Расчет собственных значений и функций Кона-Шэма, включая несколько пустых состояний.
• 3. Задается параметр optdriver=3, равный трем для расчета восприимчивости. Создается матрица для обратных значений диэлектрической проницаемости, которая сохраняется в файле с расширением SCR.

ndtset 3

# определение параметров для получения KSS файла

kssforml 3

nbandkssl 35 # число зон в файле KSS

nbandl 40

nbdbufl 5

tolwfrl 1.0d-18

# расчет экранирования

acell 3*10.217

rprim 0.0 0.5 0.5

• 0.5 0.0 0.5
• 0.5 0.5 0.0

запрещенной зоны кристаллического кремния (tgwl_l.in)

4. Задается параметр optdriver=4 для расчета собственной энергии матричных элементов, заданных набором к-точек, для получения энергии GW квазичастиц.

Пример конфигурационного файла для расчета в приближении GW представлен на рисунке 5.11. В результате расчета рассчитывается корректировка ширины запрещенной зоны (0.628 эВ), что в сумме с значением для ширины запрещенной зоны 0.594 эВ, полученной ранее, дает 1.222 эВ [1].

5.3.6 Расчет спектра поглощения с использованием уравнения

Бете-Солпитера

В данном примере рассчитывается макроскопическая диэлектрическая функция кристаллического кремния, включая экситонные эффекты, с использованием уравнения Бете-Солпитера (BSE).

На рисунке 5.12 изображено различие в расчетах оптических спектров с различными приближениями в ABINIT:

• 1) DFT в приближении LDA неправильно определяет зонную структуру и не учитывает взаимодействие электрона и дырки,
• 2) в приближении GW корректно рассчитывается зонная структура, но не

учитывается взаимодействие электрона и дырки,

Conduction bands

Valence

bands

Neutral excitations as a combination of single.

Excitonic levels

from poles of BSE

Real two-body propagator and many-body problem.

Рисунок 5.12 - Различие между подходами DFT, GW, BSE для расчета оптических спектров [1]

3) с использованием BSE корректно рассчитывается зонная структура, учитывается взаимодействие электрона и дырки.

СП

и

ф

>

•г|

optdriver 99

Рисунок 5.13 - Схема расчета с использованием BSE в программном комплексе

ABINIT [1]

Расчет в приближении BSE в программном комплексе ABINIT выполняется следующим образом (рис. 5.13):

• 1) Расчет в основном состоянии для нахождения собственных значений ?^KS и волновых функций Ч^^К8 Кона-Шэма, используется приближение LDA.
• 2) вычисление диэлектрической проницаемости выполняется после того, как получена поляризуемость для независимых частиц , используется приближение

RPA (optdriver=3). Также рассчитывается экранирование и поляризуемость для системы независимых квазичастиц, производится файл SCR.

3) выполняется стандартное приближение GW, чтобы получить волновые функции 4^/_i^QP и собственные значения энергии ^^QP для квазичастиц (optdriver=4).

ndtset 4

# определение сетки к-точек

kptopt 1 # автоматическая генерация к-точек

ngkpt 444 #

nshiftk 1

shiftk 0.0 0.0 0.0 #

# набор данных 1: расчет в основном состоянии

tolvrsl 1.0d-8

prtdenl 1

# набор данных 2: генерация KSS файла для смещеной к-сетки iscf3 -2 getden3 1 tolwfr3 1.0d-8 kssform3 3 nbandS 15

nbandkss3 10 nbdbuf3 5

chksymbreak3 0 shiftk3 0.11 0.21 0.31

Рисунок 5.14 - Конфигурационный файл для генерации WFK и SCR файлов для кристаллического кремния (tbs_l.in)

4) производится расчет с использованием уравнений Бете-Солпитера, рассчитывается макроскопическая диэлектрическая проницаемость (optdriver=99).

optdriver 99

irdwfk 1

irdscr 1

bs_calctype mbpt_sciss bs_exchange_term bs_coulomb_tenn bs_coupling bs_loband nband

bs_f reqjnesh

bs_algorithm

bs_haydock_niter bs_haydock_tol zcut

• 1
• 0.8 eV # корректировка зонной структуры

I # расчет обменного члена

II # расчет Кулоновского члена

• 0 # Приближение Тамма-Данкова
• 2 8
• 0.02 eV
• 2 # Хэйдок метод
• 200 # Максимальное число итераций
• 0.05 0 # предел сходимости
• 0.15 eV

kptopt 1 ngkpt 444 nshiftk 1 shiftk 0.11 chksymbreak 0 ecutwfn 8.0 ecuteps 2.0 inclvkb 2 icutcoul 3

0.21 0.31

acell 3*10.217

rprim 0.0 0.5 0.5

• 0.5 0.0 0.5
• 0.5 0.5 0.0

ntypat 1 znucl 14 natom 2 typat 1 1 xred

• 0.0 0.0 0.0
• 0.25 0.25 0.25 ecut 8 istwfk *1 nstep 50 diemac 12.0

Рисунок 5.15 - Конфигурационный файл для расчета макроскопической диэлектрической функции кристаллического кремния, включая экситонные эффекты, с использованием уравнения Бете-Солпитера (tbs_2.in)

Расчет выполняется в два этапа.

На первом этапе генерируются два файла WFK (волновых функций) для симметричной и смещенной k-сетки, а также файл матрицы экранирования (SCR). Конфигурационный файл представлен на рисунке 5.14.

На втором этапе рассчитывается спектр поглощения в приближении Тамма-Данкова с использованием уравнений Бете-Солпитера (рисунок 5.15). Для этого используют предварительно рассчитанные файлы WFK и SCR.

Наиболее важные результаты сохраняются в следующих выходных файлах: tbs_2o_RPA_NLF_MDF - спектр, полученный с энергиями Кона-Шэма;

tbs_2o_GW_NLF_MDF - спектр, полученный с энергиями в приближении GW;

tbs_2o_EXC_MDF - файл содержащий данные о макроскопической диэлектрической функции с экситонными эффектами, полученный с использованием уравнений Бете-Солпитера.

Рисунок 5.16 - Зависимость мнимой части диэлектрической проницаемости от частоты в различных приближениях

На рисунке 5.16 изображены спектры поглощения, полученные тремя различными подходами. Спектр RPA-KS - амплитуда первого пика занижена.

Результаты RPA-GW с корректировкой энергии не показывает значительное улучшение по сравнению RPA-KS: спектр RPA-GW просто смещается в сторону более высоких частот. Спектр BS (Бете-Солпитера) приводит к изменениям, как в оптическом пороге, так и в амплитуде первого пика. Первый пик в спектре поглощения кремния имеет сильный экситонный характер, который неправильно описан в RPA [1].

5.3.7 Расчет энергии электронного возбуждения и силы осциллятора с использованием подхода Касиды

ndtset 2

• ?DATASET 1 SCF tolwfrl 1.0d-15 nbandl 5 prtdenl 1 getwfkl 0
• ?DATASET 2 TDDFT iscf2 -1 tolwfr2 1.0d-9 nband2 12 getden2 1 getwfk2 1

acell 6 2*5 Angstrom

boxcenter 3*0.OdO diemac 1.OdO diemix 0.5d0 ecut 25 ixc 7 kptopt 0 natom 2 nbdbuf 0 nstep 25 ntypat 1 typat 1 1 xcart -0.54885 0 0 0.54885 0 0 Angstrom

znucl 7

Рисунок 5.17 - Конфигурационный файл для расчета энергии электронного возбуждения и силы осциллятора молекулы азота с использованием подхода Касиды в пределах нестационарной теории функционала плотности (ttddft_l.in)

На рисунке 5.17 представлено содержание конфигурационного файла для расчета энергии электронного возбуждения и силы осциллятора молекулы азота. Первый набор данных соответствует расчету в основном состоянии только для занятых энергетических уровней. Во втором наборе данных производится вычисления в рамках нестационарной теории функционала плотности, были добавлены еще семь незанятых уровня энергии. Коробка задается значительно большей по размеру, чем размеры молекулы азота.

В результате вычислений будет получен список энергий возбуждения, рассчитанных из разности занятых и незанятых собственных значений Кона-Шэма. Будет получены энергии возбуждения и силы осциллятора для синглетных и триплетных энергий возбуждения. По результатам можно построить зависимость поляризуемости от частоты и спектр поглощения молекулы азота [1].

5.3.8 Расчет зависимостей линейного отклика и нелинейного отклика второго порядка от частоты

В данном примере рассчитываются зависимости линейного отклика и нелинейного отклика второго порядка от частоты кристалла GaAs. Расчет выполняется в два этапа.

На первом этапе рассчитываются матричные элементы импульса -конфигурационный файл представлен на рисунке 5.18. Он имеет шесть наборов данных. В первом наборе производится стандартный самосогласованный расчет в основном состоянии, рассматриваются только занятые состояния. Во втором наборе производится несамосогласованный расчет с дополнительными незанятыми зонами, k-точки ограничены в неприводимой зоне Бриллюэна. В третьем наборе рассчитываются волновые функции для полной зоны Бриллюэна. В четвертом, пятом и шестом наборе данных вычисляются матричные элементы импульса по трем осям.

На втором этапе производится вычисление зависимостей линейного отклика и нелинейного отклика второго порядка от частоты, используя полученные данные матричных элементов импульса. Конфигурационный файл представлен на рисунке 5.19. расчет производится с помощью программы optic: optic < toptic_2.files > log &. Далее приведены входные переменные программы optic.

ndtset 6

? Набор данных 1

nbandl 4

nstepl 25 kptoptl 1 nbdbufl О

prtdenl 1 getdenl 0 getwfkl 0

?Набор данных 2

iscf2 -2

nband2 20 nstep2 25 kptopt2 1

getwfk2 1 getden2 1

?Набор данных 3

iscf3 -2

nband3 20 nstep3 25 kptopt3 3

getwfk3 2 getden3 1

?Набор данных 4

?Набор данных 5

?Набор данных 6

?Общие данные для всех наборов данных nshiftk 4

amu 69.72 74.9216

diemac 10.0 ecut 2.00 ixc 3

natom 2 nbdbuf 2 ntypat 2

rprim 0 .5 .5 .50 .5 .5 .50

xred 3*0.00d0 3*0.25d0 tnons 72*0.0

typat 1 2 tolwfr l.e-20 znucl 31 33

Рисунок 5.18 - Содержание конфигурационного файла для расчета матричных элементов импульса кристалла GaAs (toptic_l.in)

&FILES

ddkfile_l = ’toptic_lo_DS4_lWF7’, ddkfile_2 = *toptic_lo_DS5_lWF8*, ddkfile_3 = ’toptic_lo_DS6_lWF9’, wfkfile = ’toptic_lo_DS3_WFK’

/

&PARAMETERS broadening = 0.002, domega = 0.0003, maxomega = 0.3, scissor = 0.000, tolerance = 0.002

/

&COMPUTATIONS

num_lin_comp = 1, lin_comp = 11, num_nonlin_comp = 2, nonlin_comp = 123,222, num_linel_comp = 0, num_nonlin2_comp = 0,

/

Рисунок 5.19 - Содержание конфигурационного файла для расчета зависимостей линейного отклика и нелинейного отклика второго порядка от частоты (toptic_2.in)

ddkfile_l - имя файла по направлению. Подставляется предварительно рассчитанный файл матричных элементов производной по волновому вектору в данном направлении.

wfikfile - имя файла собственных значений.

domega - шаг частоты в Хартри.

maxomega - максимальная частота в Хартри.

scissor - параметр в Хартри для коррекции ширины запрещенной зоны.

tolerance - допустимое отклонение в Хартри.

num_lin_com - количество компонент из девяти линейного тензора диэлектрической проницаемости для расчета линейного отклика.

lin_comp - номера компонент линейного тензора диэлектрической проницаемости.

num_nonlin_comp - количество компонент из 27 нелинейного тензора диэлектрической проницаемости для расчета нелинейного откликаnonlin_comp - номера компонент нелинейного тензора диэлектрической проницаемости.

• 80.0
• 60.0
• (Л

Q.

Ф

• 40.0
• 20.0 -
• 0.0
• 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0

Энергия (эВ)

Рисунок 5.20 - Зависимость мнимой части линейного тензора диэлектрической проницаемости от частоты (линейный отклик)

60.0 -1

-40.0 ^J

Энергия (эВ)

Рисунок 5.21 - Зависимость реальной части линейного тензора диэлектрической проницаемости от частоты (линейный отклик)

В результате расчета создадутся выходные файлы, в которых содержатся данные о зависимостях мнимой и реальной частях линейного и нелинейного тензора диэлектрической проницаемости. На рисунках 5.20 и 5.21 приведены зависимости мнимой и реальной части линейного тензора диэлектрической проницаемости от частоты, которые получены в результате выполнения конфигурационных файлов, изображенных на рисунках 5.18 и 5.19.