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

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 10, 2014

Программный комплекс для решения задач параметрической идентификации потенциалов межатомного взаимодействия

К.К.Абгарян, М.А.Посыпкин

Аннотация— Для решения задач параметрической идентификации потенциалов межатомного

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

однокомпонентным кристаллам с ковалентным типом химической связи. Описывается многокомпонентный программный комплекс для решения задачи параметрической идентификации потенциала.

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

I. Введение

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

моделирования кристаллических структур с использованием многомасштабных подходов [1, 2]. При таком подходе, результаты первопринципного [3,4] моделирования атомной структуры и электронных свойств систем, состоящих из 200-1000 атомов используются в качестве входных данных для моделирования более сложных структур, состоящих из 1000-1000000 атомов. Такие подходы дают широкие возможности для предсказательного моделирования структур с дефектами, позволяют рассматривать различные динамические процессы (протекающие во времени), такие как диффузия, адгезия и другие. На каждом уровне масштаба можно сформулировать задачу в экстремальной постановке и применять соответствующие методы оптимизации. На первом уровне масштаба в рамках теории функционала электронной плотности [2,3] решается оптимизационная задача нахождения конфигурации базисных атомов, соответствующей минимальному значению полной энергии системы путем минимизации функционала энергии системы. В известных пакетах прикладных программ, таких как VASP

(http://cms.mpi.univie.ac.at/vasp/), Quantum ESPRESSO и других, реализован данный подход.

Работа выполнена за счет средств проекта Российского научного фонда № 14-11-00782 в ВЦ РАН.

К.К.Абгарян работает в ВЦ РАН ( e-mail: kristal83@mail.ru).

М.А.Посыпкин работает в ИППИ РАН (e-mail: mposypkin @ gmail.com).

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

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

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

II. Постановка Задачи Оптимизации Оптимизационная задача может быть формально записана как задача минимизации следующего функционала

F(f) = (f,(f)-min,fe X, (1)

,=1

в котором f, - эталонное значение i-й характеристики,

f,(f)- значение характеристики, полученное в результате расчетов для заданного набора базисных атомов, f е Rn - вектор подбираемых параметров, Щ -весовой коэффициент. Допустимое множество

X С Rn является параллелепипедом, границы которого выбираются таким образом, чтобы заведомо содержать возможный диапазон параметров. В задаче (1) требуется

найти набор параметров f е Rn, минимизирующий значение функции F(f). Такой набор будет

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

В работе предложен подход, основанный на

14

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 10, 2014

воспроизведении упругих свойств кристалла. Он

базируется на рассчитанном с помощью первопринципного моделирования значении

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

рассматриваемого материала. Количество слагаемых в (1) варьируется, в зависимости от материала.

В результате можно составить функцию вида:

F (& ) = X U(A (£>- A)2 ) ^ min (2)

i

Решение ищется на множестве X € Rn, состоящем из всех допустимых значений искомых параметров, где

£ = (£,.£ ) € Х € R - вектор подбираемых

параметров. Формирование вида минимизируемой функции F(^1,..., ^n), производится исходя из знаний о химическом составе, кристаллической структуре и типе связи моделируемого материала. Даже в простейшем кубическом случае функция F(^1,..., ^n) является многоэкстремальной.

Для решения поставленной задачи на начальном этапе необходимо определить допустимые интервалы

поиска параметров потенциала: [&, &], i= 1,...n, где n

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

потенциала Терсоффа [6]. Для нахождения точки со значениями параметров (^1,...,^n), достаточно близкими к глобальному минимуму используются методы глобальной оптимизации, такие как метод Монте-Карло, метод сканирования, Granular Random Search [6] и т.д. Далее, найденная точка используется в качестве

начальной для более точного поиска локального минимума, например методом градиентного спуска с адаптивным выбором шага, GRS (granular radial search) и другие.

Необходимо отметить, что описанные выше процедуры проводятся при определенном,

фиксированном положении базисных атомов

рассматриваемой кристаллической структуры.

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

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

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

III. Численный эксперимент. решение прикладной задачи

Рассматривается задача подбора параметров

потенциала Терсоффа[6], применительно к

однокомпонентным кристаллам. Данный потенциал является примером многочастичного потенциала,

основанного на концепции порядка связей: сила связи между двумя атомами не постоянна, а зависит от локального окружения. Такие потенциалы могут использоваться для описания свойств кристаллов с ковалентной связью (например, углерода, кремния, германия и др.).

Потенциал Терсоффа позволяет вычислить

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

E = —V V

^ coh ^ / -< у ij ’

2 i* j

Vij = fc (rij ) (VR (rij ) - bijVA (rij ) ) ,

fc(r )=< 2

1, r < R - R„

f

1 - sin

{n(r - R) ^

V 2Rcut JJ

R - Rcut < r < R + Rcut,

0, r > R + RcUt,

V r )

S^-"J eXP (-e2S (rij - re ))

VR (rij ) =~Т^ eXP (-в2^ (rj - re ))

Va (rtJ)

S -1 SDe S -1

r -r ) ,

'ij e

f

exp

V

-eS2(r« - re )

bj =(1+(*„ )n Г'2 n.

Z = X fc r) g (j)»«,

k* i, j

®ik = eXP(^3 (rij - rik )3 ),

g J)=1+fd) - d2+(hc2 . )2

V d J d + (h - cosj)

Переменные rj означают расстояния между атомами

с номерами i и j соответственно, через J обозначен

угол, образованный векторами, соединяющими атом i с атомами j и к соответственно. Потенциал Терсоффа включает 12 параметров, специфичных для моделируемого

15

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 10, 2014

вещества: De, re,в, S, n,y,A, c, d, h, R, Rcut. При этом, параметры R и Rcut определяются исходя из

полученных экспериментально геометрических характеристик вещества и не нуждаются в подборе.

Далее рассматривается задача идентификации параметров потенциала Терсоффа для изотропных кристаллов с кубической кристаллической решеткой. Согласно закону Гука, малые деформации пропорциональны напряжениям [7].

а = Се, (3)

где а - тензор напряжений, е - тензор деформаций, С -тензор упругих констант, компоненты которого определяются выражениями

С-

1 д2 Е

V де tде: ’

1 j

i,j=1,...6

(4)

где е- - соответствующие компоненты тензора деформации кристалла.

Таким образом, С- (С) - значения упругих констант,

рассчитанных по формуле (4) с использованием значения когезионной энергии, вычисленной с помощью выбранного потенциала межатомного

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

упругих констант. Число независимых упругих констант зависит от кристаллической структуры моделируемого материала. Т ак, в случае, если моделируется кубический кристалл остается лишь три независимые упругие константы С11, С12, С44 . Оставшиеся компоненты

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

f (С) = Ж (Е(С) - Е)2 + (a(£) - a)2 +

+ щ(в(С) - B)2 + ^4(сЧС) - С')2 + ,

+ Ж (С44С) - С44 )2 + ^6 (с (С) -Z )2,

где C = (De, re, в, S, n,y,A, c, d, h) - вектор

идентифицируемых параметров. В целевую функцию входят следующие величины

Е - когезионная энергия системы, отнесенная к числу атомов моделируемой совокупности (удельная энергия);

a - константа решетки (характерный размер);

B - модуль объемной упругости;

С' - модуль сдвига;

С44 - постоянная эластичности;

Z - параметр Клейнмана, который характеризует

дополнительное изменение решетки при сдвиговой деформации. Это изменение является следствием стремления атомов занять энергетически оптимальное положение.

Величина С44 является компонентом тензора

упругих постоянных, а модули объемной упругости и сдвига вычисляются через компоненты С11 и С12 по формулам

B = (Сц + 2С12)/3, С' = (Сц - С12)/2.

Эталонные значение компонент целевой функции, полученные экспериментально либо

квантомеханическими методами для монокристалла кремния, взяты из работы [2]:

Е = -4.63, a = 5.43, B = 0.97,

С = 0.51, С44 = 0.79, zf = 0.52.

Параллелепипед

X = [СС = Rn : С — xi — С} для определения

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

С = (0.5, 0.5, 0.5, 0.5, 0.1,5 ■ 10-8, 0.5, 10000, 1, - 2) С = (10, 5, 5, 5, 2, 3 ■ 10-6,3, 200000, 30, -0.1).

Все расчеты производятся для одной ячейки кристаллической решетки кремния, состоящей из 18 атомов (Рис. 1). При этом учитываются взаимодействия не только между атомами ячейки, но и с соседними атомами в пределах, определяемых радиусом отсечения

Rcut ‘

Рис 1. Одна ячейка кристаллической решетки кремния

Когезионная энергия рассчитывалась в соответствии с формулами вычисления потенциала Терсоффа, после чего выполнялось деление вычисленного значения на число атомов в ячейке для получения удельной энергии. Компоненты B, С' вычислялись методами

аппроксимации отношения изменения энергии к величине соответствующих малых деформаций решетки. Характер деформации, применяемой для вычисления постоянной эластичности С44 таков, что

часть атомов смещается с позиций, определяемых деформацией, занимая энергетически-оптимальное

16

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 10, 2014

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

Был реализован программный комплекс,

включающий в себя следующие модули:

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

2. Модуль вычисления основных характеристик материала и целевой функции.

3. Модуль оптимизации, выполняющий поиск минимума на допустимом множестве.

Рассмотрим алгоритмы и подходы, составляющие первый модуль и лежащие в основе заполнения координат кристаллической решетки [8]. При описании кристаллической структуры рассматриваемого материала, соответствующей заданной химической формуле первоначально задаются свойства симметрии (с помощью одной из 230 федоровских групп симметрии) и позиции Уайкова для базисных атомов. Задание федоровской группы симметрии позволяет построить элементарная ячейку рассматриваемого материала, с помощью групп симметрии, ее можно транслировать (размножить). Элементарная ячейка

задается системой неравенств, а группы симметрии -системой уравнений. Задание позиций Уайкова для базисных атомов (разрешенных, в рамках

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

Благодаря симметриям кристалла для отыскания оптимальной структуры (соответствующей минимуму внутренней энергии системы) достаточно рассмотреть конечное множество всех пар атомов {(i,j)} таких, что i пробегает множество неэквивалентных позиций базисных атомов в пределах одной элементарной ячейки, а j- множество позиций атомов, в пределах этой и соседней с ней ячейках.

Пусть с помощью системы неравенств задана элементарная ячейка соответствующая федоровской группе симметрии (a, b, c - параметры, задающие элементарную ячейку). Для заданной химической формулы в рамках заданной группы можно построить множество различных конфигураций базисных атомов, находящихся в соответствующих позициях Уайкова. Рассмотрим k-ю конфигурацию базисных атомов

(к = 1, ~ ) :

Xk = {Х1 (Х1, ?1, Z1X...X N (xN, yN, zN )}eW ) Здесь N - число базисных атомов, W- множество всех возможных конфигураций. Для каждого базисного атома из заданной конфигурации необходимо определить разрешенные, в рамках федоровской группы симметрии позиции Уайкова. В общем виде, для каждого базисного атома, координаты центра могут

--к ---к _

быть заданы X i (xt, yt, zt) = X (i). Обозначим через

0 --к ------0

W множество всех X (i)GW i . Ограничения

накладываемые принадлежностью к конкретной позиции Уайкова имеют следующий вид:

X = C}a, yt = C?b, Zi = C3c (5)

где 0 < С] < 1, 0 < Сi < 1, 0 < С3 < 1

Для заданной федоровской группы симметрии существует конечное число операций симметрии

(j = 1, K , K - число операций симметрии данной федоровской группы), которые позволяют размножить базисные атомы внутри элементарной ячейки, а затем построить образы атомов соседних с ней ячейках. Обозначим через

X (i, j) = {x(i, j), y (i, j), z(i, j), i = 1, N,j = 1, K}

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

Xk (i, j) = St (X1 (i)) + Т~}

Здесь Si - матрица коэффициентов, Т j -вектор.

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

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

(i = 1N,j = 1K), {X(i, j )}gW

Для каждого атома-образа соответствующего атома из W0 при таком построении будут выполняться

аналогичные (5) ограничения, но с учетом применения операций симметрии.

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

W = UWi,(i = 1N,j = 1K).

i, j

Программный модуль с помощью которого происходит «размножение» базисных атомов внутри элементарной ячейки, а затем их трансляция в соседние ячейки с использует информацию по федоровским группам симметрии, соответствующим операциям симметрии и позициям Уайкова , которая

17

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 10, 2014

хранится в специально созданной базе данных.

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

Модуль вычисления основных характеристик материала достаточно прост и основан на формулах потенциала Терсоффа. Наиболее сложным компонентом программного комплекса является модуль оптимизации, реализованный в рамках библиотеки BNB-Solver[9]. Решаемая задача относится к классу «черного ящика», т.к. явное выражение для целевой функции не известно и ее значение есть результат работы сложной вычислительной процедуры. Для решения данной задачи был применен подход, состоящий в сочетании методов глобализации поиска, обеспечивающих диверсификацию поиска в пределах допустимого параллелепипеда и локальных методов, позволяющих находить локальный минимум в заданной окрестности.

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

Были опробованы два локальных оптимизационных алгоритма. Первый алгоритм - метод градиентного спуска с адаптивным выбором шага. Использовалась численная аппроксимация градиента методом конечных разностей. Второй алгоритм GRS (granular radial search) предложен в [10]. Метод работает следующим образом: из начальной точки выполняется случайный сдвиг по одному из параметров в заданном диапазоне g . Если удается уменьшить значение целевой функции, то вновь найденное значение замещает начальное и выполняется новый сдвиг. Начальное значение при этом выбирается случайно в параллелепипеде X . Если процент «удачных» сдвигов, приведших к уменьшению значения целевой функции, становится меньше заданного порога, то уменьшается гранулярность g. Алгоритм останавливается или по превышению верхней границы количества итераций или когда гранулярность становится слишком малой.

Результаты предварительных вычислительных экспериментов показали, что алгоритм градиентного спуска работает хуже для данной задачи по сравнению с GRS. При ограничении на число итераций и точность, дающих приблизительно одинаковое время работы, алгоритм GRS обеспечивает существенно более высокую точность. В эксперименте с 8-ю запусками из различных начальных приближений с помощью метода GRS было найдено оптимальное значение 0.1987, тогда как применение градиентного спуска дало минимальное значение целевой функции, равное 0.0002 .

Вычислительный эксперимент был организован следующим образом. Были сгенерированы 1024 случайных начальных приближений, равномерно

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

четырехядерными процессорами Intel(R) Xeon(R) CPU E5620 @ 2.40GHz. С целью ускорения вычислений

запускалось 8 параллельных процессов, каждый из которых обработал 128 начальных приближения. Общее время расчетов составило примерно 4 часа.

Из найденных локальных минимумов дальнейшей обработке подвергались только те, значение целевой функции в которых не превосходило 0.001 . Для таких точек производилась дополнительная численная

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

U(х) = {у е Rn :| х{ — yi \<р|, Р = 0.1. Количество точек вычисления целевой функции было положено равным 1000.

Для 13 сгенерированных точек, удовлетворяющих условию f (4) ^ 0.001 проверка на локальный

минимум прошла успешно. Был проведен анализ отклонения найденных параметров от эталонных значений. По каждому из 10 параметров было вычислено отклонение di от эталонного значения по

формуле:

d =

\ х— r \

\ r \

i = 1,...,10,

где - значение параметра, r - эталонное значение

параметра. Точки на Рис. 2 соответствуют значениям отклонения для 10 определяемых параметров.

Dt гг Р S п у Я с d h Рис. 2. Отклонения в значениях 10-ти параметров для найденных локальных минимумов

Полученные результаты показывают, что с заданной целевой функцией удается достаточно точно определить

параметры De, re ,fi, n, h . Остальные параметры

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

18

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 10, 2014

IV. Применение методов параллельных и

РАСПРЕДЕЛЕННЫХ ВЫЧИСЛЕНИЙ

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

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

Также технологии параллельных и распределенных вычислений можно применить для исследования точности определения параметров. Для этого необходимо провести большое число экспериментов и построить распределение отклонения параметров по большой выборке. Поскольку вычисления абсолютно независимы, для этой цели идеально подойдут распределенные вычисления [11]. В частности, перспективным представляется организация так называемых добровольных вычислений на платформе BOINC[12]. В этой модели вычислений узлами распределенной системы являются компьютеры пользователей, которые подключаются по сети Интернет к серверу проекта и получает от него задания для вычислений.

V. Заключение

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

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

необходимо провести вторичный фитинг,

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

Также планируется расширение списка применяемых оптимизационных алгоритмов, распараллеливание

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

Библиография

[1] Абгарян К.К. Применение оптимизационных методов для моделирования многослойных полупроводниковых наносистем // Труды Института системного анализа Российской академии наук. Динамика неоднородных систем. 2010 г. Т.53(3). С.6-9.

[2] K. Abgaryan. Multiscale computational model of multy-layer semiconductor nanostructures. /IV International Conference on Optimization Methods and Applications. (OPTIMA- 2013). Petrovac, Montenegro. September 2013.Abstracts. pp.13-14.

[3] P. Hohenberg, W. Kohn, Phys. Rev. 136, B864, 1964.

[4] W. Kohn, L. J. Sham, Phys. Rev. 140, A1133, 1965.

[5] K.Abgaryan, D.Bazanov, I.Mutigillin. The simulation of the nitridation of the sapphire surface (001). Materials XXI International Conference "Interaction of ions with the surface"VIP 2013, Т.2. С. 99-103.

[6] J.Tersoff, Empirical interatomic potential for silicon with improved elastic properties. Phys. Rev. B. 1988, V.38, P.9902-9905.

[7] Ч.Киттель. Введение в физику твердого тела. Из-во техникотеоретической литературы. Москва.1957 г.

[8] Абгарян К.К., Хачатуров В.Р., Компьютерное моделирование устойчивых структур кристаллических материалов. Ж. вычисл.матем. и матем. физ., 2009. Т. 49. № 8. C. 1517-1530.

[9] Посыпкин М.А. Методы и распределенная программная инфраструктура для численного решения задачи поиска молекулярных кластеров с минимальной энергией //Вестник нижегородского университета им. Н.И. Лобачевского № 1. 2010. С. 210-219.

[10] Powell D. Elasticity, Lattice Dynamics and Parameterisation Techniques for the Tersoff Potential Applied to Elemental and Type III-V Semiconductors : дис. - University of Sheffield, 2006.

[11] Афанасьев А. П. и др. Развитие концепции распределенных вычислительных сред //Проблемы вычислений в распределенной среде: организация вычислений в глобальных сетях. Труды ИСА РАН.-М.: РОХОС. - 2004. - С. 6-105.

[12] Anderson D. P. Boinc: A system for public-resource computing and storage //Grid Computing, 2004. Proceedings. Fifth IEEE/ACM International Workshop on. - IEEE, 2004. - С. 4-10.

19

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 10, 2014

Software for solving problems of parametric identification of the interatomic interaction

potential

K.K. Abgaryan, M.A. Posypkin

Abstract - Various optimization methods are applied for parametric identification of the interatomic interaction potential. We consider the problem of selecting the parameters of the Tersoff potential, applied to single-component crystals with covalent type of chemical bond. The paper describes multi-component software package for solving the problem of parametric identification of the potential.

Keywords - optimization methods, parametric identification, the objective function, gradient methods, software package

20