Современная компьютерная реализация прямого расчета постоянной Маделунга кристаллов типа ab

УДК 53.072

И.Е. Еремин, М.С. Сычев

СОВРЕМЕННАЯ КОМПЬЮТЕРНАЯ РЕАЛИЗАЦИЯ ПРЯМОГО РАСЧЕТА ПОСТОЯННОЙ МАДЕЛУНГА КРИСТАЛЛОВ ТИПА AB

Рассматривается эффективный алгоритм расчета постоянной Маделунга, разработанный на базе метода векторно-матричного описания структуры кристаллической решетки в совокупности с методом Харрисона. Проводится анализ результатов вычислительного эксперимента для кристаллов типа AB.

Постоянная Маделунга, метод Харрисона, прямой метод расчета, дипольный момент, условная сходимость.

I.E. Eremin, M.S. Sychev

RECENT COMPUTER REALIZATION OF THE DIRECT CALCULATION MADELUNG

CONSTANT OF AB-TYPE CRYSTALS

The effective algorithm of calculation Madelung constant is considered.

This algorithm is based on the method of vector-matrix presentation structure of crystal lattice and Harrison method. The results of computing experiment for AB-type crystals are analyzed.

Madelung constant, Harrison method, direct method of calculation, dipole moment, conditional convergent.

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

В физике конденсированного состояния большое внимание уделяется величине постоянной Маделунга - AM . Она характеризует силовое взаимодействие между элементами кристаллической решетки, следовательно, и ее стабильность. Знание данной константы позволяет определить полную энергию кристаллической решетки, а также модуль упругости. В работе [1] был представлен метод векторно-матричного описания структуры кристаллической решетки и предложен модифицированный алгоритм прямого расчета постоянной Маде-лунга. Проведенный вычислительный эксперимент на базе разработанной компьютерной программы, использующей данный алгоритм, показал, что его применение возможно только для расчета постоянной Маделунга кристаллических решеток обладающих скомпенсирован-

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

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

Совсем недавно У. Харрисоном был предложен оригинальный способ расчета постоянной Маделунга [2], не страдающий от условной сходимости, следовательно, гарантирующий получение правильного результата. Ключевая идея метода Харрисона состоит в следующем: исходная элементарная ячейка кристалла разбираемого типа одновременно расширяется в трех измерениях, т.е. применяется метод расширяющихся кубов, после чего вводится дополнительное граничное условие - сфера радиуса R, вписанная в расширяющийся куб и исключающая электростатическое взаимодействие F ионов, оказавшихся снаружи, обусловленное значением соответствующих решеточных сумм. Очевидно, что суммирование значений F ионов, оставшихся внутри сферы Харрисона, приводит к колебаниям результирующего заряда (который может быть либо положительным, либо отрицательным), величина которого нарастает с увеличением радиуса R. Поэтому, для компенсации неустойчивости суммарного заряда внутри сферы, к рассмотрению первоначально полученного результата добавляется учет тонкой оболочки с радиусом R и зарядом Q, электростатический вклад которой считается равным - Q/R , где Q - результирующий заряд, полученный при суммировании зарядов всех ионов, находящихся внутри сферы Харрисона, включая заряд центрального иона.

Таким образом, способом численного определения постоянной Маделунга, наиболее подходящей для автоматизации ее расчета может считаться метод Харрисона. Однако его существующие программные реализации оказываются достаточно громоздкими, поэтому исследователи, работающие в затрагиваемой области знаний, поступают следующим образом. Внутренняя энергия решеток, обладающих скомпенсированным суммарным зарядом, рассматривается в рамках прямого расчета AM , а более сложные случаи - с помощью метода Харрисона.

С целью реализации возможности создания универсального вычислительного средства, предназначенного для расчетов постоянной Маделунга кристаллических решеток произвольных структур, авторами был разработан следующий алгоритм их автоматизации. Отметим, что в его основу были заложены ранее предложенный способ компактного описания кристаллической решетки [1], а также общая идея метода Харрисона.

Во-первых, задаются элементарные матрицы-генераторы нечетных Zнечет и четных

Zчет координационных слоев, отражающие расположение заряженных частиц в решетке конкретного типа. Для примера будут рассмотрены структуры МаС1 и СвС1, для них матрицы-генераторы имеют следующий вид:

Z

Z

та

-1 +1 Z^NaCl +1 -1

+1 -1 ’ чет -1 +1

cscl

0 0 ZCsCl _ +1 0

0 -1 , чет 0 0

Во-вторых, вводится значение Ь номера внешнего (предельного) координационного слоя, на основании которого формируются двухмерные векторно-матричные модели ZЬ и

ZL-1 расположения зарядов исследуемой трехмерной кристаллической структуры.

В-третьих, на основании введенного значения Ь в автоматическом режиме формируются матрицы, эффективно характеризующие количество однотипных позиций пространственных узлов внешнего КЬ и предшествующего ему КЬ-1 координационных слоев вида, обладающие размерами Ь х Ь и (Ь - 1)х(Ь -1).

В-четвертых, посредством поэлементного умножения матриц ZL на КЬ и ZL-1 на КЬ-1 рассчитываются глобальные матрицы, соответственно, внешнего МЬ и предшествующего ему МЬ-1 координационных слоев.

В-пятых, учитывая перевод условной линейной величины межъядерного расстояния конкретной кристаллической структуры к значению длины ребра а элементарного куба, определяется граничное условие метода Харрисона - радиус R ограничительной тонкой оболочки, с помощью формулы:

Я _ Ь  к ,

где к - коэффициент упомянутого перевода линейных величин, равный для решетки типа МаС1 единице, а для решетки типа С8С1 составляющий 3-0,5.

В-шестых, выполняется вычисление решеточных сумм для каждого из заданных координационных слоев, начиная с первого, матрицы которых формируются из глобальных матриц вида МЬ или МЬ-1, путем удаления необходимого количества их внутренних строк и столбцов. При этом применяется граничное условие метода Харрисона, проверяющее соотношение

кд/

при выполнении которого подобающие элементы решетки учувствуют в расчете значения текущего члена ее ряда Маделунга, проводимого по формуле:

л _ £ £

1_1 ._1 V У

В-седьмых, каждая из полученных решеточных сумм сохраняется в памяти ЭВМ для участия в определении итоговой величины постоянной Маделунга.

В-восьмых, посредством суммирования величин зарядов, находящихся внутри сферы радиуса Я, включая заряд центрального иона, вычисляется величина компенсирующего заряда Q.

В-девятых, на базе сохраненных данных рассчитывается конечное значение постоянной Маделунга по формуле:

Ь

Ам _ £ Л - .

I _1

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

Для проведения общей оценки эффективности предлагаемого алгоритма был разработан пробный вариант программы расчета постоянной Маделунга, реализованный на макроязыке вычислительной среды Ма1:ЬаЬ, выбор которой обусловливался использованием векторно-матричных моделей, а также высоким быстродействием проводимых вычислений.

С помощью названного программного продукта, был проведен вычислительный эксперимент. Расчеты были выполнены для трех разновидностей кристаллических структур: гранецентрированной кубической решетки типа МаС1 (рис. 1а), объемно-центрированной кубической решетки структурного типа С8С1 (рис. 1б) и алмазоподобной кристаллической структуры а^пБ (рис. 1в). Результаты вычислительного эксперимента расположены в таблице. Для выявления точности результатов разбираемых вычислений, а также оценки време-

ни их производительности tFa6 было рассмотрено большое количество координационных

слоев L, включающих в себя весьма значительные числа частиц N .

Отметим, что для проведения разбираемого эксперимента использовался компьютер, оснащенный процессором Intel Core 2 Duo CPU E4500, тактовая частота ядер которого составляет 2,20 и 2,21 ГГц, а объем ОЗУ - 3,24 Гб.

Рис. 1. Элементарные ячейки NaCl, CsCl и a-ZnS

Анализ результатов вычислительного эксперимента, позволяет сделать вывод, что применение выше изложенного алгоритма автоматизированного расчета постоянной Маделунга, даёт возможность рассчитывать значение константы, для кристаллических решеток, не обладающих скомпенсированным зарядом. Кроме того, полученные результаты, с одной стороны, соответствуют общепринятым значениям, с другой стороны, обладают высокой точностью. Также необходимо отметить, что использование метода векторно-матричного описания кристаллической решетки сокращает массив непосредственно обрабатываемых данных почти в 200 раз, что приводит к существенному увеличению скорости производимых вычислений.

Результаты вычислительного эксперимента

К ристаллическая решетка Примечание

структурный тип №С1 структурный тип CsCl структурный тип о^пБ

Результаты вычисли- тельного экспери- мента Ам = 1,74822296405540 ^ = 1 мин. 0,1 сек. Ам = 1,76344636372008 ^ = 1 мин. 0,4 сек. Ам = 1.63752776623190 ^ =20,65 сек. L = 500 N = 1 003 003 000

Ам = 1,74793679225170 ^ = 8 мин 7,1 сек. Ам = 1,76190306809904 = 8 мин 11,4 сек. Ам = 1.63847926830246 ^ =2 мин. 42,3 сек. L = 1000 N = 8 012 006 000

Ам = 1,74790238753186 ^ = 26 мин. 16,2 сек. Ам = 1,76255923407786 = 26 мин. 11,2 сек. Ам = 1.63803451177794 ^ =9 мин. 1,2 сек. L = 1500 N = 27 027 009 000

Ам = 1,74764245654470 ^ = 1 ч. 02 мин. 13,8 сек. Ам = 1,76296658864009 ^в = 1 ч. 02 мин. 14,9 сек. Ам = 1.63828454255004 ^ =22 мин. 42,3 сек. L = 2000 N = 64 048 012 000

Ам = 1,74758061114039 ^ = 2 ч. 01 мин. 36,8 сек. Ам = 1,76271905735291 ^ = 2 ч. 01 мин. 22,3 сек. Ам = 1.63795149128957 ^ =41 мин. 38,2 сек. L = 2500 N = 125 075 015 000

Справочное значение AM = 1,7476 Ам = 1,7627 Ам = 1,6381 Лит. источник [3]

Основываясь на выше изложенном можно резюмировать, что предлагаемый авторами алгоритм автоматизации расчетов постоянной Маделунга, обладает весьма неплохими характеристиками. Он позволяет рассчитывать значение Ам не только для кристаллических решеток, имеющих нулевой дипольный момент, но и обладающих не скомпенсированным зарядом. Также его практическая реализация в виде компьютерной программы показала, что результаты расчетов имеют высокую точность, а время, затраченное на их получение не вели-86

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

ЛИТЕРАТУРА

1. Еремин И.Е. Модифицированный алгоритм прямого расчета постоянной Маделунга / И.Е.Еремин, М.С.Сычев // Информатика и системы управления. 2010. № 3(25). С. 27-34.

2. Harrison W.A. Simple calculation of Madelung constant / W.A. Harrison // Physical Review B. 2006. V.73. P. 212103.

3. Киттель Ч. Введение в физику твердого тела / Ч. Киттель. М.: Наука, 1978.

Еремин Илья Евгеньевич -

кандидат физико-математических наук, доцент кафедры «Информационные и управляющие системы» Амурского государственного университета

Сычев Михаил Сергеевич -

аспирант кафедры «Информационные и управляющие системы» Амурского государственного университета

Статья поступила в редакцию 23.07.11, принята к опубликованию 14.11.11