Учет дальнодействия при распределенном молекулярно-динамическом моделировании систем большой размерности Текст научной статьи по специальности «Физика»

Учет дальнодействия при распределенном молекулярно-

динамическом моделировании систем размерности

Воронова Л.И. fvoronova2001@mail.ru), Тен Э.А.

Курганский государственный университет Работа выполнена при поддержке РФФИ, проект № 01-07-96506

Введение

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

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

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

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

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

За последнее десятилетие рядом авторов [1-4] были разработаны подходы, используемые в локальных программных комплексах [5-11], позволяющие моделировать конденсированные системы содержащие ~103 частиц. Однако в настоящее время актуальна задача разработки вычислительных моделей систем большой размерности (число частиц 104- 105 и более) с использованием новых информационных технологий.

В Курганском госуниверситете на кафедре прикладной математики и компьютерного моделирования под руководством д.ф.-м.н, профессора Л.И.Вороновой разработана информационно-исследовательская система (ИИС) «Шлаковые расплавы». Работа проводилась при поддержке Российского Фонда Фундаментальных Исследований по направлению «Создание и развитие информационных, вычислительных и телекоммуникационных ресурсов для проведения фундаментальных исследований» [12].

ИИС «Шлаковые расплавы» это совокупность информации, физико-химических и математических методов и моделей, технических, программных, других технологических средств и специалистов, предназначенная для генерирования специализированной научной информации, ее обработки и принятия на ее основе прогнозных решений.

ИИС построена на базе новых информационных технологий ^ORBA, WEB, XML) как распределенная система, с удаленным доступом через Интернет, интегрирующая доступные

информационно-вычислительные и телекоммуникационные ресурсы (набор ПК и низкоскоростные каналы связи).

ИИС обеспечивает реализацию компьютерных экспериментов (КЭ) в рамках комплексной модели многокомпонентного шлакового расплава для систем большой размерности (104-106 частиц), хранение, автоматизацию обработки и исследования модельных результатов. [13].

Методика учета дальнодействия

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

Фундаментальной проблемой МД-метода, в применении к сильновзаимодействующим системамявляется адекватное описание потенциала межчастичного взаимодействия.

Приближение для потенциала определяет вид математической модели. В случае ионной модели основой является ион с характеристиками: т1 - масса, q1 - заряд, о1 - радиус жесткой сферы, Г -радиус-вектор, у1 - скорость; р(г^) - парный, сферически симметричный потенциал,

сила,

и (г1 ,г2 ,--,гп) = ^р(гд) - потенциальная энергия системы, =^-[дрдг ) ~

действующая на частицу.

Существует много разновидностей модельных ионных потенциалов, применяемых при МД-моделировании, воспроизводящих основные особенности потенциальной кривой [14-18].

В потенциале р^ц-) обычно выделяют два основных вклада. На расстояниях, меньших г0

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

дальнодействия: р(г^ ) = ркор (г^ ) + ркул(г^ ).

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

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

Кулоновская энергия системы N частиц в модельном кубе с длиной ребра Ь и его бесконечных репликах рассчитывается по формуле

1 N N q.q ,

и^ , (1)

2 •_А • _А Г

п 1 = 1 ] = 1 Ц,п

где п =( П], п2 , п3) - вектор трансляций модельного куба, П], п2 , п3 - целые числа.

Расстояние между частицей в модельном кубе и частицей в реплике, определяется как:

г ~ = Ь',п

а)

б)

Рис. 1. Трансляции модельного куба Рис. 2. Распределения плотностей зарядов

Суммирование ведется от исходной ячейки п = (0; 0; 0) по всем ячейкам, окружающих

базовую, до бесконечности. При этом слагаемое при i = ] опускается, когда п = 0.

Ряд (1) медленно сходится, поэтому в ИИС потенциальная энергия рассчитывается по методике Эвальда [19], первоначально разработанной для решеточного суммирования кулоновского взаимодействия в ионных кристаллах. Для этого на электростатическое поле (рис. 2а), создаваемое

дискретными зарядами частиц в базовой ячейке с плотностью р^(г) = qi5(r -г{) накладывается

искусственное поле (рис. 2б), которое создавалось бы этими частицами, если бы они имели гауссово распределение зарядов противоположного знака:

р (;) = qiа3 ехр(-а2 (г - гг)2)

Полный заряд иона в рассматриваемой точке гг = 0 в этом случае равен:

а ехр\ -а г \ да да да з/2/222

\Р,Сгт = |-^^& = | ¿х | аУ { а ехр(-а (х_ +у +' ,&=

V

V

-да -да

= qi-^ [ ехр (-а2 х2 )) (ах)-^ [ ехр (-а2 у2 )) (ау)-^ [ ехр (-а2 ¿2 ) (аz) = qi, л/^ уп

* -да * -да * -да

да

где | ехр(-2 )& = 4п - интеграл Эйлера-Пуассона.

-да

Полный потенциал в данной точке решетки будем искать в виде суммы

кул пр , обр

р =щ +

(2)

двух различных, но связанных между собой потенциалов рРР и р(бр . Смысл разделения потенциала

пр обр

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

пр обр

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

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

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

2

-да

Vj =

rr

перекрываются с областью рассматриваемой точки (рис. 2а). Вклад в этот потенциал от каждого иона решетки 1 имеет вид

ь - ±.¡р ()У- | рМ у. (3)

г1 Г] У. У\У£ т

Первое слагаемое в (3) обусловлено точечным зарядом, второе — зарядом с гауссовым распределением, заключенным внутри сферы £ радиуса т, описанной вокруг узла решетки 1, третье — распределением заряда вне указанной сферы.

Рассмотрим вычисление интегралов в формуле (3): 1. У. - объем, ограниченный сферой £ радиуса т^

3 I 2 "" 2 \ а ехр\-а т )

¡Ру (т)йУ = qj | ---йУ = {переходя к сферическим координатам, получим:}

V.£ У.£

rj -3 Рsina Ъ. 24j-g-rj -а■ + 2q■«rj

га exp \-а р i 2 , Г ■ п jn Г j ч -a j ¿-Hi г -t2

= q ■ I-_ и L p dp I sinQ dQ I dp =--J ._ J e 1 I e dt.

0 Vn3 0 0 4n 0

2. V\VS - все пространство без объема VS:

р1 (r) dV = q I a exp{-a2 r 2 )dVГ = q ■ exp (-a2 pP ) pdP sinQ d6\dm = Iqa e~«2 r1 2

I 1>dV = qj I a exp^a r 'dV = qJ \a exp\~3 P /pdpI sinQdQ |dp

V\VS r V\VS yn ■ r r ■ Vn о 0

л[ж

Подставив полученные выражения в формулу (3), получим:

{ . ат у Л ( ат, Л

qj i

Vj =—— rr 1 1

2qj -a2 r12 + 2 qj f;e-t2

4П 4П

{ e-t2 dt -ЩО-e-«2 rj2 = q-L ! 2 / e-'2 dt

-'«•■ J Vn r -I™- J

0

Vn

0

x

= —ет/с(ату), где ет/с(х) = 1 - ет/(х) = 1 —,= Iе и йи -дополненная до единицы функция ошибки. Таким образом, потенциал в прямом пространстве имеет вид:

ррр = Ър} = X —ет/с (ату), где 1 * 1. 1=1 1=1 т1

Выражение для потенциала в прямом пространстве можно записать с помощью вектора трансляций модельного куба п =( П], п2 , п3):

N ет/с(а ■ т - )

рПР =1^1-^ .

- т -

п 1 -1 1] ,п

Потенциал частиц роб]Р, создающих искусственно накладываемое поле, соответствует

непрерывному гауссовому распределению (рис. 2б), поэтому его несложно вычислить в обратном пространстве с помощью преобразования Фурье.

Разложим р0р и плотность зарядар в ряды Фурье:

0ф =2Cleikr, p=%p^ekr

k к

где к = 2п ■ ^Ь'т'Ь= Ь — вектор обратной решетки, И = ( 1,т,р ), 1,т, Для нахождения коэффициента Ск используем уравнение Пуассона:

V2робр =-4пр . (4)

Подставляя робр и р в формулу (4), получим:

_2 -- --

iк г а _ х-* ik г к

с-гк ■е = 4п2 Рке ^ ст = ■

_ к ^к ~ — ~к -»2 • к к к

Полная плотность распределения зарядов в модельном кубе, где распределение обусловлено содержащимися в ней узлами и «хвостами» распределений всех узлов пространства вне модельного куба (рис. 3 а) равна полной плотности зарядов всего пространства, содержащего только узлы модельного куба (рис. 3б), т.е.

_ егк7 = 2 qjа3 ехР(-а2 (г - ГJ) 2) (5)

к 2 гт ()

j=1 \п

а) б)

Рис. 3. Распределение плотности заряда в модельном кубе

Умножим обе части равенства (5) на е 1кг и проинтегрируем по объемам:

рк{ е>кг ■е-кгаЛ = \ 2

Л V j=1

• 2

N 3 -а2 (г - г )2

а.а е \ 1 ' .т-

"1 „-1 к г

е-ikг ¿V ^рг Л = е 4а ^ а

N

■ е

—i ■ к ■ г

1=1

к

В = — ■е к Л

4а

N

•2 а-'

■=1

-Vк■ г ■

^ „4а2 N --

4п е ^ -i■ к■ г/

ск А -у ' 2 1е 1

к

Л

к __

с 4а2 N . к ~ -— 4п е 4а2 N . к ( )

2^2—2 "1-екг1 ■е = — е--2 а.-е1'кГ-г1)

обр 4п х-* е р = — ■-2-^ 24

к к 1 =1

к к 1 =1

В рассматриваемой точке i потенциал будет иметь вид:

обр Х-* 4 п е р =2—л

4а2 N ■2

Ь(г, -г1)

к ~ к 1=1

Мы не должны учитывать вклад в потенциал, обусловленный распределением заряда в самой рассматриваемой точке i (рис. 2б пунктирная линия).

,3 _ I „2 ^ 2 '

{в^м = а ехр(~а2 г 2) ¿V=^.

л г 13

п3 ■ г

V V

Таким образом, потенциал в обратном пространстве имеет вид:

к2

4а2 N

обр х-* 4п е i■ к ■ (г1 - г 1 ) 2 "

р. =2 А '2 —=

а

4П

к - л 1=1

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

2

к

1

2

к

U кул = U

кул(пр)

+ U

кул(обр) т j кул(пост)

U

кУл ( пР ) 1

N

1

U

NN

=2 2 q (7 = - 22 qi 2qj

erfc(a■ ~ )

ijn

1 N 1

Uкул (обр) = 2 2 q. ■ (pfp = 1 2

2 i=1

4n

exp

' i=1 ( k2 ^

i=1 j=1

iJn

v a у

Nq 2 qj icos [[ - rj)]+i-sin [[ (ri- rj )])=

2 k V

i=1 j=1

(6) (7)

= 1 2 4n e

4a2 N

N

2 qi 2 qj

2^ V k2 ^ ^ ^ -»j

2 k V k i = 1 j =1

■cos

k- (ri- rj)] где qi 2 qj ■sin

i=1 j=1

k■ (r-rj )]= 0-

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

/

ехр

( k2 ^

1 N N

Uкул (обр) = 1 22 qi 2 qj л

2 k i = 1 j = 1 ^

4n

v 4a^ у

■ cos

■(r - rj )],

1 N

U кул (пост) = 2 q..

2 f=1

1 N 2 2 N пост = 1 у 2 qi a = a y q2

(i = о 2 г- = г 2 qi •

2 i=1 Vn Vn i=1

(8) (9)

Уравнения Эвальда (7) - (9) могут быть упрощены для ускорения вычислений. Значение

параметра а можно выбрать так, чтобы в прямом пространстве учитывались только члены с п = 0, то есть суммирование вести по модельному кубу, и увеличить число членов в сумме по обратному пространству.

При п = 0 получаем:

/ ^ ^ егк(а-г, )

икул (пр) = 2 чг2 ч,——- (10)

I I <] гу

Двойное суммирование по г и по , в обратном пространстве (8), ведущее к квадратичной зависимости №2 объема вычислений от размера моделируемой системы, последовательно приводится к двум суммам порядка N следующим образом:

ехр

( к2 ^

U кул (обр) = 2п Х-*

U = V2 k2

k * 0 k

v 4a у

2

2 cj + 2 sj

vj у v j у

(11)

где cj = qj cos(k-rj ), Sj = qj sin(k- rj ).

Сила, действующая на частицу может быть получена из уравнений (10) и (11), дифференцированием по координатам. После преобразований получим

-~кул -~кул (пр) -¿кул (обр)

Fi = Fi + Fi

где

F

кул (пр)

N

=qi 2 qj

i * j

2a exP (-a2 ri2) + erfc(a ■ rij)

кул ( обр)

F i =— 2 k ■

i V

k *0

vn

exp

r

iJ

ri

ij

( k2 ^

v 4ay

2 cj ■si - 2 s,

v J J v j у

c

(12)

(13)

2n-

при k =-h формула (13) примет вид:

L

r

2

2

2

r

r

F

кул (обр) 2

exp

I h

п2 h 2 ^ а2 L2

I/ ■ si - I

V j Vj )

Т2 -- Ь 2

^ Ь * 0 Ь

Слагаемое цк^л(пост в уравнении (6) является константой и следовательно не участвует в суммировании сил, однако вносит существенный вклад в общую кулоновскую энергию.

Обсуждение результатов

Для существенного увеличения размеров модельных систем и производительности вычислений в ИИС «Шлаковые расплавы» применяется параллельная реализация МД-алгоритма, с распределением расчетов на основе СОЯБА -технологии. Распределенные данные представляются в виде набора вычислителей - однотипных объектов выполняющих расчет потенциалов и сил одновременно на всех ЭВМ в вычислительной сети. При этом возможно масштабирование вычислительной мощности ИИС, причем, время расчета больших систем сокращается пропорционально количеству задействованных ПК. Оптимизация параметров моделирования осуществляется автоматически.

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

Уравнения (9) - (14) зависят от двух параметров: а (параметр относительной сходимости рядов)

- 2п - 3

и ктах =-Ьтах (максимальный вектор обратной решетки). В системах с N > 10 расчет в прямом

Т

пространстве ~ N и является наиболее затратным по времени. Для уменьшения вычислительной нагрузки вводим параметр гкул - радиус обрезания, обеспечивающий суммирование по прямому пространству только в сфере с этим радиусом.

Минимизация гкул возможна при увеличении а. Как следует из уравнений (10), (12) при а ^ да ег/е(а-х) ^ 0. С другой стороны, для убыстрения сходимости по обратному пространству а нужно

(

уменьшать, т.к. из (11), (14) следует, что при а ^ 0 ехр

.2 Л

4а'

^ 0.

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

Рис. 4. Зависимость кулоновских сил при фиксированном параметре а от:

2

а) г в прямом пространстве; б) х = Ь в обратном пространстве.

2а ехр(-а2г2) ег/е(а - г ) --- +-~-

Модуль силы в прямом пространстве при а >0 Fпр(а,r) =—¡=--

Ып r

является

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

r

min

max

(F пр (a,r))

= s.

(15)

Для убывающей с расстоянием функции min(fпр (а,r))= Fпр {(х,rкул ),

max(Fпр(а,r))= Fпр(а,rmin), где rmin - минимальное расстояние между частицами (сумма радиусов жестких сфер).

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

длины вектора h (рис. 4б).

exp

( 2 ^ п x

F (а, x) = ^

v а2L2 j

L

, где x = h .

(16)

*2

=s.

(17)

Задавая параметры а и е можно найти ктах .

обр {а, х)) = Р обр [а, ^ ах) тах( обр (а, х)) = ""Р°^(ад) Для вычисления сил с точностью не хуже 1% можно выбрать е =10-3.

В табл.1 приведены значения параметров процедуры Эвальда системы ЫаС1 для различных гкул

при фиксированном е, в зависимости от числа частиц в базовой ячейке.

Параметры процедуры Эвальда системы №С1

Табл.1

x

N 103 203 303 403 503

L, Ä 28 56 84 112 140

гкул = 8 Ä, а = 3,14-109

h2 max 28 80 139 201 262

^куфр) c 0,2 3 12 24 52

^кул(обр) c 0,3 22 188 895 3301

t, c 0,5 25 200 919 3353

гкул = 10 Ä, а = 2,36-109

hi2 max 18 53 95 140 187

tKул(пр) c 0,2 8 22 48 91

c 0,3 16 107 445 1413

t, c 0,5 24 129 493 1504

гкул = 12 Ä, а = 1,88-109

h2 max 13 37 68 103 139

c 0,4 15 40 97 204

tкул(обр) c 0,3 10 80 345 874

t, c 0,7 25 120 442 1078

гкул = 14 Ä, а = 1,55-109

hi2 max 9 28 52 78 107

tKул(пр) c 0,4 29 64 156 303

!КУл(о6р) c 0,3 6 52 181 666

t, c 0,7 35 116 337 969

гкул = 16 Ä, а = 1,31 -109

h2 max 7 21 40 61 85

}кул(пР) с 0,5 42 112 248 575

^куфЬр) с 0,2 4 34 138 493

г, с 0.7 46 146 386 1068

гкул = 18 А, а = 1,13 • 109

И2 тах 6 17 32 49 68

^кул(пр) с 0,6 45 208 368 848

^кул(°бр) с 0,2 3 20 130 351

г, с 0.8 48 228 498 1199

гкул = 20 А, а = 0,98 -109

И2 ах 5 14 26 40 56

^кул(пр) с 0,8 46 217 591 936

^кул(°бр) с 0,2 2 21 141 254

г, с 1 48 238 732 1190

Здесь: гкул"рр, гкул°рр _ время расчета кулоновской энергии и силы по прямому и обратному пространству за один шаг, г _ время моделирования одного шага.

Тестирование проводилось в локальном варианте на РепйишШ 500МГц.

Анализ полученных данных, приведенных в табл.1 показывает, что в небольших системах (К~103) для гкул = 8 А и гкул = 10 А время счета (~0.5 с) минимально и примерно поровну делится между прямым и обратным пространством. Дальнейшее увеличение радиуса обрезания приводит к ухудшению результата.

Для систем содержащих ~8-103 частиц значения радиуса обрезания можно выбирать в пределах от 8 до 12 А. При этом временные затраты по прямому и обратному пространству «перераспределяются», однако общее время счета ~ 25 с. Дальнейшее увеличение гкул значительно увеличивает общее время счета.

Для больших систем, содержащих ~ 10 частиц и более оптимальное значение для ткул = 14 А. Уменьшение радиуса обрезания - существенно увеличивает время счета по обратному пространству, при увеличении радиуса обрезания - резко возрастает время счета по прямому пространству.

Например, в системе с Ы=503 при выборе небольшого радиуса обрезания (гкул = 8 А) требуется в несколько десятков (~60) раз больше времени для расчета сил и энергии по обратному пространству, чем по прямому, а при гкул = 20 А = Ь/7 время моделирования по прямому пространству почти в 4 раза превышает время по обратному.

При увеличении радиуса обрезания, время счета по прямому пространству монотонно возрастает, по обратному - монотонно убывает. Их сумма дает общее время моделирования, которое достигает оптимума при определенном гкул. Для систем, состоящих из 103 - 2.7-104 частиц оптимальным будет гкул = 10А, для систем, состоящих из 104 _-105 частиц - гкул =14 А.

Выводы

Для реализации молекулярно_динамического моделирования сильновзаимодействующих систем большой размерности в ИИС «Шлаковые расплавы» применяется распараллелеливание алгоритма с помощью технологии СОЯБА.

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

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

Проведены тестовые эксперименты по определению оптимальных значений параметров на системе №С1 в зависимости от числа частиц в системе и радиуса обрезания в прямом пространстве. Из результатов экспериментов следует, что:

_ для систем малых размеров (К<103) в силу периодических граничных условий радиус обрезания следует выбирать порядка половины ребра расчетного куба гкул = Ь/2;

_ для больших систем (N>10*) радиус обрезания следует выбирать гкул ~ 14 А. Такое значение радиуса обрезания существенно ограничивает число расчетов парных взаимодействий в прямом

пространстве и ускоряет сходимость сумм. При этом обеспечивается точность расчета сил ~ 10-3, что вполне приемлемо для МД-моделирования;

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

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

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

Список литературы

1. Alder B. J., Wainwright T. E. Studies in Molecular Dynamics. I. General Method // J. Chem. Phys., 1959. V.31. №2. P.459-466.

2. Rahman A., Stillinger F. H. Molecular dynamics study of liquid water // J. Chem. Phys., 1971. V. 55. №7. P.3336-3359.

3. Лагарьков А.Н., Сергеев В. М. Метод молекулярной динамики в статистической физике // УФН, 1978. Т.125. Вып.3. С.409-448.

4. Beeman, D. Some multistep methods for use in molecular dynamics calculations // J. Comput. Phys., 1976. V.20. P.130-139.

5. Eastwood J.W., Hockney R.W., Lawrence D.N. P3M3DP - the three-dimensional periodic particle-particle, particle-mesh program // Comp.Phys.Comm., 1980. V.19. №.2. P.215-261.

6. Soules T. F. A molecular dynamics calculation of the structure of sodium silicate glass // J. Chem. Phys., 1979. V.71. №11. P.4570-4578.

7. Mitra S. K., Hockney R. W. Microheterogeneity in simulated soda silica glass // The structure of noncrystalline materials. New York, 1982. P.316-325.

8. Белащенко Д.К. Компьютерное моделирование структуры и свойств некристаллических оксидов // Успехи химии, 1997. Т.66. №9. С.811-844.

9. Зильберман О.П., Черников А.Н., Зильберман П.Ф., Знаменский В.С. Исследование контактного плавления в системе KCL-NaI-KI // Материалы 5-го Российского семинара "Компьютерное моделирование расплавов и стекол". Курган, 2000. С.107-108.

10. Бойко Г.Г., Паркачев В.А. Изучение структуры фосфатных стекол методом молекулярной динамики // Тр. 2 Российского семинара "Компьютерное моделирование физико-химических свойств стекол и расплавов". Курган, 1994. С.6-9.

11. Voronov V.I., Voronova L.I., Rizhov N.I., Gusev A.I. Properties of a System FeO-SiO2-CaO-MgO by Results of Computer Experiment // Proceedings Second International Conference on Mathematical Modeling and Computer Simulation of Metal Technologies, "MMT-2002". College of Judea and Samaria, Israel, 2002. P.1-76 - 1-83.

12. Воронова Л.И., Рыжов Н.А, Воронов В.И., Тен Э.А., и др. Подсистема распределенного молекулярно-динамического моделирования информационно-исследовательской системы "Шлаковые расплавы". Свидетельство об отраслевой регистрации разработки № 3158. Зарегистрировано в Отраслевом фонде алгоритмов и программ 04.02.2004. Государственный координационный центр информационных технологий, МО РФ.

13. Воронова Л.И., Воронов В.И., Рыжов Н.А., Гусев А.И., Тен Э.А. Информационно-исследовательская система "SLAG MELT // Тр. IV Международной конференции "Компьютерное моделирование 2003". Санкт-Петербург, 2003г. С.264-267.

14. Mulliken R.S. Magic formula, structure of bond energies and isovalent hybridization // J.Phys.Chem., 1952. V.56. №3. P.295-317.

15. Fumi, F. G., and Tosi M. P.. Ionic Sizes and born repulsion parameters in the NaCl type alkali halides // J. Phys. Chem. Solids, 1964. V.25. P.31-43.

16. Woodcock L.V., Angell K.A., Cheeseman P. Molecular dynamics studies of the vitreous state: simple ionic system and silica // J. Chem. Phys., 1976. V.65. №4. P.1565-1577.

17. Angell C.A., Clarke J.H.R., Woodcock L.V. Interaction potentials and glass formation: a survey of computer experiments // Adv.Chem.Phys., 1981. V.48. P.397-453.

18. Воронова Л.И., Бухтояров О.И., Вяткин Г.П. Расчет параметров потенциала Me-O методом MNDO для МД моделирования в ионно-ковалентном приближении: 1. Анализ применимос ти MNDO-расчетов для построения потенциальных кривых // Расплавы, 1994. №6. С.50-57.

19. Sangster M.J.L., Dixon M. Interionic potentials in alkali halides and their use simulation of molten salts // Adv.Phys., 1976. V.25. №3. P.247-342.