Определение электропроводности жидкого цезия в широком диапазоне температур методом рекурсии с использованием параметризации Андерсена Текст научной статьи по специальности «Физика»

УДК 539.2

ОПРЕДЕЛЕНИЕ ЭЛЕКТРОПРОВОДНОСТИ ЖИДКОГО ЦЕЗИЯ В ШИРОКОМ ДИАПАЗОНЕ ТЕМПЕРАТУР МЕТОДОМ РЕКУРСИИ С ИСПОЛЬЗОВАНИЕМ ПАРАМЕТРИЗАЦИИ АНДЕРСЕНА

А.Н. Соболев, АЛ. Мирзоев

В работе предложено усовершенствование методики расчета электропроводности жидких металлов и сплавов методом рекурсии. Расчет выполнен для кластера жидкого цезия для 9 различных температур. Проведено сравнение полученных результатов с экспериментом.

Введение

В кристаллических металлах вычисление электронной структуры упрощается из-за периодичности атомной структуры (то есть периодичности атомных потенциалов). В случае же веществ, не обладающих трансляционной симметрией, к которым относятся неупорядоченные сплавы, стекла, а также расплавы и растворы, расчеты электронной структуры с использованием зонной теории становятся практически невыполнимыми. С другой стороны, успешное развитие химии твердого тела, основанное на представлении о локальном характере межатомных взаимодействий, показывает, что свойства вещества в основном определяются лишь взаимодействием атомов с соседями нескольких ближайших координационных сфер. В связи с этим в последнее время для расчетов электронной структуры неупорядоченных веществ исследователи все чаще обращаются к использованию метода сильной связи, основное достоинство которого состоит в локальности описания и применимости к системам с любым типом беспорядка, включая жидкие металлы и сплавы. Трудности метода заключаются в необходимости иметь небольшой, но полный базисный ряд функций, которые были бы одновременно весьма локализованными. Обычно при расчете матрицы гамильтониана Ни оставляют только двухцентровые члены и трактуют эти интегралы как подгоночные параметры, в силу чего метод долгое время рассматривался как полуэмпирический. Однако в последнее время взгляд на возможности метода был пересмотрен. В большой степени это связано с возможностью построения удобного набора базисных функций на основе идей, развитых в методе линейных маффин-тин орбиталей (JIMTO) O.K. Андерсона [1].

Представление гамильтониана в методе сильной связи приводит к очень разреженной матрице HSJ (большое количество матричных элементов равны нулю), поскольку, если приближение сильной связи справедливо, то орбиталь на данном узле взаимодействует только с ближайшими соседями. Это обстоятельство позволяет сформулировать быстрый и удобный алгоритм диагона-лизации и нахождения локальной и полной плотностей электронных состояний неупорядоченных систем, развитый Р. Хэйдоком и получивший название метода рекурсии [2]. Метод рекурсии оказался одним из наиболее удачных методов моделирования неупорядоченных материалов и широко используется для расчета электронных и фононных спектров жидких металлов [3,4]. Настоящая работа посвящена применению метода рекурсии к расчету электропроводности жидкого цезия в широком диапазоне температур по формуле Кубо - Гринвуда [5,6]. Основная проблема, возникающая при этом, - сложность точного определения собственного вектора гамильтониана для собственного значения, равного энергии Ферми данной системы. Решение данной проблемы было предложено Крамером и Вейром в 1978 году [7]. Эти авторы предложили для нахождения собственного вектора использовать процедуру фильтрации.

Однако сходимость метода для фильтрующего оператора, предложенного в работе [7], очень слабая (для достаточной точности необходимо, по крайней мере, 3-105 итераций), поэтому такая процедура требует большого количества машинных ресурсов. В данной работе предлагается усовершенствование метода фильтрации, существенно снижающее объем вычислительных затрат.

В первой части работы дается формальное описание используемых методов, во второй части приводятся полученные результаты и их обсуждение.

Соболев А.Н., Мирзоев А.А.

Определение электропроводности жидкого цезия е широком диапазоне температур методом рекурсии...

Используемые методы

Метод рекурсии. Метод рекурсии позволяет численно определить диагональный элемент функции Грина С0(Е) и затем вычислить локальную плотность состояний по формуле

п(Е) = —- Im G0 (Е + is). к

(i)

Основная идея метода состоит в том, что для заданного в виде матрицы «-го порядка гамильтониана Н строится ортогональный базис, в котором он имеет тридиагональную форму в виде матрицы Якоби:

Н -

а,

О

а-

О Ъ;

О О

h а3

v

(2)

Для тридиагонализации гамильтониана используется следующая процедура. Вначале в качестве исходного вектора берется орбиталь на выбранном узле в базисе сильной связи, для которого нужно вычислить плотность состояний. Если обозначить через г/о^-мерный вектор начального приближения, то остальные векторы находятся по рекуррентным формулам:

Ь0

и0\и0

<*п =

< ип \Н\ип >

и

п+1

(.H-an)iin-bnun_u

<un\un>

h2 _ < un+1 I un+1 >

л+i - ;

<un\un>

(3)

Здесь коэффициенты ап и Ъп имеют смысл вероятностей: ап - вероятность для электрона остаться на месте при действии гамильтониана, Ъп - вероятность перейти на какую-либо из других орби-талей. В нашей работе используется приближение сильной связи. Это значит, что электрон в результате действия гамильтониана может перейти только на атомы из ближайшего окружения. Действие второй и последующих сфер соседей пренебрежимо мало.

Исходный вектор щ состоит из нулей, кроме одного элемента, равного 1, который соответствует выбранному узлу. Тогда

От{Е)=<щ (Е-НТ1 щ >= ¿о КЕ~ао~Ь\ /

(4)

!{Е - ах - Ъ\ /[Е ~ а2-----!{Е - - Ь^Е))•••]})

где Ооо(Е) - нулевой диагональный элемент функции Грина, а функция ¿(Е) называется терминатором и служит для обрывания бесконечной цепной дроби на определенном уровне, который называется уровнем рекурсии.

Расчет электросопротивления. Согласно формуле Кубо - Гринвуда [5, 6], диагональные элементы тензора проводимости при нулевой температуре в представлении собственных функций выглядят следующим образом:

e2h 2

ajj = ТГ X Н Ет Iъ I Еп >1 SiEm - EF )siEn ~Ef)> а тм

(5)

где Cía - объем образца, h - постоянная Планка, энергия Ферми и у,- у-й компонент оператора скорости.

Представив ¿-функцию в виде мнимой части функции Грина и использовав выражение

^3(Е-Ет)<Ет\/(Е)\Ет >=g(E){f(E)}Em=E ,

(6)

где {}\Е означает среднее по собственным векторам с энергией Е, g(E) - плотность состояний в энергии Еу уравнение (5) может быть записано в следующей физически прозрачной форме:

Vjj =—n(EF)D(EF).

(7)

Здесь С1а - объем атома, п(ЕР) - плотность состояний на атом, В{ЕГ) - функция перескока

ОДр) - -А Иш 1т{< Ет \v.GiEp + ¡ф, Ет >}Ет=

£-> О

(8)

В{Е)<) может быть численно вычислена как проекция средней локальной плотности состояний на состояния и мы можем посчитать этот множитель, используя метод рекурсии. Собственные векторы \Ет> вычисляются с использованием техники фильтрования, описанной ниже. Оператор скорости вычисляется по соотношению:

Получение собственных векторов |у> с достаточной точностью методом фильтрации - часть вычисления, требующая наибольшего времени.

Процедура фильтрации

В работе [2] в качестве фильтрующего оператора предлагался оператор

(Н~а)(Н-Ь)

(10)

где а и Ъ - точки, в интервале между которыми находится А — а —Ь — величина интервала между а и Ъ\ Еь - энергия Ферми.

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

(х-аХх-Ь) (11)

Д2

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

Как уже отмечалось, сходимость такого метода очень медленна (при Л = а-Ь = 20 эВ для достаточной точности необходимо по крайней мере 3-105 итераций), такая процедура требует большого количества машинных ресурсов. В данной работе предлагается метод нахождения собственного вектора с заданной энергией, не требующий вычислительных затрат. Он основан на непостоянности точек а и 6, а также введении так называемого «машинного нуля» - а Ф 0, ограничивающего снизу квадрат каждой компоненты вектора состояния.

Алгоритм нашего метода достаточно прост. Предлагается критерий сходимости:

8 ~ \(Н — Ерегш - квадрат модуля вектора отклонения состояния от собственного на энергии

Ферми. Очевидно, что при у/, равному искомому собственному вектору, 8=0. Также очевидно, что эффективность фильтрации будет наиболее высокой, если ветви параболы фильтрующего оператора (рис. 1) пересекают ось энергии как можно ближе к границам зоны возможных значений энергий системы (то есть как можно ближе к дну зоны проводимости и энергии Ферми). Этого можно добиться, если вместо постоянства величин а и Ь, как было в [2], корректировать их положение соответствующим образом. Для этого на каждой второй итерации производится проверка на сходимость по критерию 8 < где 8о есть 8\ рассчитанная за две итерации до этого. Для выбора наиболее оптимальных параметров сужения выполнены вычисления с различной скоростью сужения; были получены результаты, которые демонстрируют, что наилучшая сходимость достигается при симметричном движении точек а и Ъ к энергии Ферми на 0,1 эВ при ¿> < ¿ь, и движении точек а и Ь от энергии Ферми на 0,3 эВ при 8> <5Ь. Кроме того, в последнем случае вектор состояния принимается равным вектору, записанному за две

Рис. 1. Вид оператора фильтрации в пространстве энергий

Соболев А.Н., Мирзоев А.А.

Определение электропроводности жидкого цезия в широком диапазоне температур методом рекурсии...

итерации до этого. Таким образом, А в операторе фильтрации по крайней мере на 0,2 эВ больше, чем разница между наибольшим и наименьшим собственными значениями разложения вектора состояния по собственным функциям гамильтониана, что обеспечивает хорошую сходимость процедуры фильтрации.

Для уменьшения машинного времени, требуемого для фильтрации, и упрощения самих вычислений вводится еще один параметр - с, так называемый «машинный нуль». При достижении квадратом компоненты вектора состояния значения е обнуляем эту компоненту. Очевидно, что чем больше е, тем меньше необходимо шагов для полной фильтрации. Но, обнулив хотя бы одну компоненту вектора состояния с ненулевым числом ближайших соседей, мы уже не можем с помощью данной процедуры получить точный собственный вектор, но только вектор, близкий к собственному, поэтому данный параметр принимается достаточно большим (а = 10*4) только на ранних этапах фильтрации. Кроме того, вычисления показывают, что значение 3= 10~3 достаточно для получения довольно точных результатов.

Полученные результаты

В данной работе проведен расчет электропроводности расплавов цезия для температур 323, 573, 773, 973, 1073, 1173, 1373, 1673 К. Список координат атомов структурных моделей жидкого цезия (1000 частиц) в широком интервале температур был получен построением структурных моделей обратным методом Монте-Карло (ОМК) [8] и был любезно предоставлен нам проф. Д.К. Белащенко. Для построения матрицы гамильтониана в методе сильной связи использовалась предложенная в [9] экспоненциальная интерполяционная формула для интегралов перескока между орбиталями соседних узлов как функция расстояния между ними. Для определения параметров интерполяционной формулы использовались интегралы перескока для ОЦК- и ГЦК-решеток цезия, полученные методом ЛМТО в работе [10]. Была проведена аппроксимация этих значений функцией

яу

Усф^^+А-е* , (12)

где - интеграл перескока между орбиталями а и р атомов / и у, зависящий от расстояния

между этими атомами; ¥®р, А, 7?о - параметры аппроксиманты; Яш - радиус Вигнера - Зейтса для расчетного кластера. Значения параметров приведены в таблице.

Таблица

Параметры аппроксимации точек, полученных Андерсеном, экспонентой

Параметр Значение Абс. погрешность

V" I 0 0

А 122,88933 ±128,35254

Яо -4,67957 ±0,57253

Нами проведен расчет локальной плотности Б-состояний, усредненной по 20 атомам кластера, и электропроводности для каждой из 9 температурных точек при различных коэффициентах параметризации гамильтониана. Для этого была построена сетка из различных пар параметров, и для каждой пары параметров рассчитывалась электропроводность. Собственный вектор для расчета последней находили путем 1000 итераций действия оператора фильтрации на произвольное состояние. Результаты расчета электропроводности представлены на рис. 3.

0 035 0,030 0.025 0 020 Ч > 001 ОЛЮ Ч 0 005 0 И

1.9

-т— 1 в

параметры Андерсена

аппроксимаций экспонентом

—г*

20

—г~

2.2

'А':

—г~ 2,4

2.6

Рис. 2- Интегралы перескока по Андерсену

3 000 000-

2 2500000 S

о

 2000000 -

«3 1 500 ООО -о

с& &

С t 000 000 â ti

g 500000-

0,

эксперимент [и] ■ теоретические данные

i—н—■—i—•—i—■—i—•—i—«—i—»—i—-—i—»—i—3 200 400 600 600 1000 1200 1400 1600 J 800 20ÖC

Температура, К

Рис. 3. Зависимость электропроводности от температуры

Расчет проводился при указанных выше условиях: параметрах сужения +0,1/-0,3 эВ, экспоненциальной параметризации гамильтониана. Результаты записывались при достижении S значения 6- 10~3. Результаты, полученные для двух последних точек, 1673 и 1923 К (они не показаны на рисунке) некорректны с физической точки зрения для всех рассчитанных параметров гамильтониана, что объясняется неприменимостью метода для очень высоких температур.

Выводы

Из рис. 3 видно, что полученная теоретическая зависимость вполне соответствует экспериментально полученной кривой, повторяет ее качественно и не намного отходит от нее в количественном отношении. Отклонения от экспериментальной кривой, по-видимому, можно объяснить просто неточностью метода либо ограничениями, связанными с применимостью параметризации функции гамильтониана экспонентой, но величина погрешности метода окупается его быстротой.

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

Литература

1. Tank R.W., Arcangeli С. An Introduction to the Third-Generation LMTO Method // Phys. stat. sol. 2000, 217, p. 89.

2. Haydock R. The Recursive Solution of the Schrödinger Equation // Solid St. Phys., 1980,35, p. 215.

3. Дунаевский C.M., Розова M.H., Кленкова H.A. Электронная структура графитовых нанот-рубок // Физика твердого тела - 1997. - т. 39. - № 6. - С. 1118.

4. Мирзоев A.A., Гельчинский Б.Р., Бушуев A.C. Некоторые возможности метода рекурсии для изучения электронных и фононных спектров конденсированных сред // Высокотемпературные расплавы. - 1996. - № 1. - С. 17.

5. Kubo R., J. Phys. Soc. Jpn., 1957,12, p. 570.

6. Greenwood D.A., Proc. Phys. Soc. London, 1958, 71, p. 585.

7. Kramer В. and Weaire D., J. Phys. C: Solid State Phys., 1978,11, p. L5.

8. Белащенко Д.М. Самосогласованность формы парной корреляционной функции плотных неупорядоченных систем // Металлы, - 1988. -№ 4. С. 172-174.

9. Anderson O.K., Pawlowska Z., Jepsen О. Illustration of linear muffin tin orbital tight-binding representation: compact orbitals and charge density in Si // Phys.Rev. B, 1986, 34, № 8. - P. 5253.

10. Andersen O.K., Jepsen O. et al. in Highlights of Condensed Matter Theory, ed. by F. Bassani, F. Fumi, M.P. Tosi.-North-Holland, NY, 1985.

11. Каган Д.Н. и др. Жидкометаллические теплоносители тепловых труб и энергетических установок. - М.: Наука, 1988. - 268 с.

Поступила в редакцию 26 апреля 2003 года