CC BY
43
ВЫЧИСЛИТЕЛЬНАЯ СЛОЖНОСТЬ АЛГОРИТМА ЭВАЛЬДА В МЕТОДЕ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
© Тен Э.А.*
Курганский государственный университет, г. Курган
В статье предлагается методика оценки вычислительной сложности алгоритма Эвальда учета дальнодействующих вкладов в межчас-тичные взаимодействия, которая позволяет существенно снизить временные затраты при молекулярно-динамическом моделировании систем, содержащих более 10000 частиц в модельном кубе.
Одной из актуальнных задач научного направления «Создание новых материалов с заданными свойствами» является развитие математического моделирования строения вещества, с использованием ЭВМ, в основе которого лежит метод молекулярной динамики, заключающийся в расчете траекторий движения частиц модельной системы (например, оксидного расплава) и позволяющий определить целый комплекс свойств (структурные, термодинамические, транспортные).
При молекулярно-динамическом моделировании приходится учитывать ряд специфических особенностей, увеличивающих время эксперимента: для изучения транспортных свойств необходимо максимально увеличить количество шагов моделирования; большая кривизна функций потенциала около минимума требует очень малого шага интегрирования уравнений движения для сохранения устойчивости. Однако основной причиной резкого увеличения времени моделирования является учет дальнодействующего вклада, пропорционального г"1, в межчастичных взаимодействиях при исследовании систем, содержащих заряженные частицы. Кулоновская энергия определяется по формуле:
1 т т т NN г-1 ,
и — % I Е Е Е Е-
где п — (п1, «2, «3) - вектор решетки с целочисленными компонентами,
г - — \гг.] + ^'п - расстояние между одной частицей в модельном кубе с г] ,п 1 1
длиной ребра Ь и другой частицей, находящейся в отображении куба (рис. 1). Суммирование ведется от исходной ячейки п — (0; 0; 0) по всем,
* Доцент кафедры «Прикладная математика и компьютерное моделирование», кандидат технических наук
окружающих ее, копиям до бесконечности (слагаемое при i = j опускает-
На текущий момент существует несколько алгоритмов, позволяющих вычислять кулоновскую энергию с различной вычислительной сложностью и точностью. Широко используемым является алгоритм Эвальда.
Метод суммирования по Эвальду был предложен впервые в [4] как методика для вычисления дальнодействующего взаимодействия между частицей и ее бесконечными периодическими образами. На электростатическое поле (прямое пространство), создаваемое дискретными зарядами N
частиц ячейки с плотностью p (r) = qt8(r - rt) накладывается искусственное поле (обратное пространство), которое создавалось бы этими частицами, если бы они имели гауссово распределение зарядов p (r) = q ai exp (-a2 (r - r)2 )/V? (a - параметр ширины гауссовой кривой) противоположного знака. В ходе вычисления суммы потенциалов по прямому и обратному пространствам вклады по зарядам с гауссовым распределением компенсируются и то обстоятельство, что распределение имеет гауссов характер исключается. Кулоновскую энергию можно найти с помощью двух быстро сходящихся рядов:
ся, когда n = 0).
где U ^ул(пр - энергия, вычисляемая в прямом пространстве; цкул(обр) - в обратном пространстве.
2я .
m m m
exp
U ^ (°бр) = ^ llm у у у _ V m^“ h^mh^mh^m
2я
4a2
( N V ( N ^
Ус + Ул
V M ) V j== )
N
a2 ■77 У q
V я ¡=\
где k = — • А Ф 0 , h = (h, h, h) - вектор с целочисленными компонен-
тами, cj = qj cos (k • rj ), s;. = qj sin
i(k • rj ).
Для вычисления сумм с заданной точностью используются параметры: гкул и |йтах | - радиусы обрезания, определяющие количество слагаемых
в суммах соответственно в прямом и обратном пространствах. Минимизация гк'л возможна при увеличении а, т.к. в прямом пространстве при а ^ ж вг/с(а-х) ^ 0. С другой стороны, для убыстрения сходимости по обрат-
х2
ному пространству а нужно уменьшать, т.к. при а ^ 0 exp I -
4a2
■0.
Для определения вычислительной сложности алгоритма Эвальда необходимо оценить количество операций для всех его фаз:
1. Вычисление потенциальной энергии и силы по прямому пространству.
2. Вычисление потенциальной энергии и силы по обратному пространству.
Пусть С], С2 - количество операций, необходимое для выполнения каждой из перечисленных фаз, тогда общее количество операций:
С = С, + С2
В прямом пространстве, при использовании радиуса обрезания, учитываются взаимодействия только частиц, находящихся на расстояниях, не больших гкул. Куб с длиной ребра Ь, при этом, разбивается на ячейки с
длиной ребра гкул, количество которых
L
. При равномерном распре-
делении частиц в модельном кубе, количество частиц в одной ячейке: N • I —— I . Так как взаимодействие частиц рассматривается только в бли-
жайших соседних ячейках, тогда:
(
С = С • N •
27 • N •
зЛ
2
k
2
2
k
r
r
L
L = :
• VN
р
С, = с, • 271 — \ • N Для обратного пространства количество векторов обратной решетки
равно (2 lAmax |), поэтому С2 = с2 • 8 |йп
• N - с3 • N.
Здесь С], с2, с3 - константы, определяемые как количество атомарных операций (сложение, умножение), требующихся для выполнения соответствующей операции.
Таким образом, слагаемые оценки количества операций (вычислительная сложность алгоритма) зависят только от количества частиц и от параметров алгоритма. Величины параметров метода Эвальда, хотя и гораздо меньше количества частиц системы, но зависят от этого количества. Например, в работе А.М. Баженова [1] выбираются параметры Г’”' = L/2,
-*■2
hmax = 5, позволяющие рассчитывать кулоновские силы в системах до 100 частиц с точностью до 1%. Для систем до 1000 частиц выбирают обычно
—2
фикированные параметры, в работе [5]: = L/2, hmax = 6.
Однако, выбор параметра гкул = L/2 приводит в прямом пространстве к расчетам ~ N2 , так как:
С = с, • N^27• N• 1 j = ^С • N2
и поэтому, в больших системах такой выбор является неприемлемым. В пакете AMBER [3] (разработан исследовательской группой под руководством П.Колмана в Калифорнийском университете США), предназначенного для МД-моделирования биомолекул, значения гкул выбирают в пределах 8-12 А. Например, результаты молекулярно-динамического моделирования системы SiO2-CaO на одном шаге для разного количества
о — 2
частиц с параметрами гкул = 8А, hmax = 25 приведены в табл. 1, на рис. 2.
Таблица 1
Время моделирования одного шага для постоянных параметров алгоритма Эвальда
N 103 5103 104 2-104 5-104 6404 7404 9-104 105
t, c 0,2 8 22 77 380 556 774 1191 1324
N
m
i=1
з
На рис. 2 просматривается квадратичная зависимость, поэтому для фиксированных параметров алгоритма Эвальда вычислительная сложность с увеличением числа частиц будет расти пропорционально квадрату.
В системах с N > 103 в работе [2] параметры находят из уравнений:
Рпр (а, гкул)
-----т------г = е
Рпр (а,гШ1П)
Р °бр (а, ¿1* )
Робр (а,1)
= 6
где ^(а, г) -модуль силы в прямом пространстве;
^(а x) -модуль силы в обратном пространстве; е - относительный порядок отброшенных членов (задается исследователем);
гты - минимальное расстояние между частицами (сумма радиусов жестких сфер).
Рис. 2. Квадратичная зависимость времени моделирования одного шага от числа частиц при постоянных параметрах алгоритма Эвальда
Кроме относительного порядка отброшенных членов е, исследователем задается радиус обрезания прямого пространства гкул, который следует выбирать в зависимости от числа частиц в молельном кубе так, чтобы время моделирования одного шага было оптимальным (равенство времени расчета сил и энергий прямого и обратного пространств). Для е = 10-4 в табл. 2, на рис. 3 приведены значения оптимальных (по скорости вычислений межчастичных взаимодействий на одном шаге) параметров алго-
ритма Эвальда и время моделирования одного шага системы 8Ю2-СаО для различных гкул, в зависимости от числа частиц в базовой ячейке.
Таблица 2
Параметры алгоритма Эвальда и время моделирования одного шага
N 103 5-103 104 2-104 5-104 6-104 7-104 9-104 105
гкул, А 8 10 10 12 12 12 14 14 14
h2 h max 26 40 58 57 95 105 81 93 99
t, c 0,2 6 14 46 172 222 270 370 395
На рис. 3 прослеживается линейная (близкая к линейной) зависимость, поэтому вычислительная сложность алгоритма Эвальда будет линейно зависеть от числа частиц только при выборе оптимальных параметров этого алгоритма.
Рис. 3. Линейная зависимость времени моделирования одного шага от числа частиц при вычисленных параметрах алгоритма Эвальда
Однако нужно еще учесть, что выбор параметров алгоритма Эвальда
является компромиссом между точностью и скоростью вычислений, так
как увеличение точности приведет к существенному увеличению времени
счета, из-за необходимости выбора больших значений радиусов обрезания
гкул и h ' -*1 ^lmax•
Список литературы:
1. Баженов А.М. Структура, теплофизические и транспортные свойства ионных расплавов и неидеальной плазмы. - Диссертация на соискание ученой степени кандидата физико-математических наук. - Свердловск: УПИ, 1983. - 110 с.
2. Воронова Л.И., Тен Э.А. Учет дальнодействия при распределенном молекулярно-динамическом моделировании систем большой размерности [Электронный ресурс] // Электронный журнал «Исследовано в России». -2004. - Режим доступа: http://zhumal.ape.relam.ru/articles/2004/208.pdf.
3. Case D.A. et al. Amber 5. Technical report. - University of California, 1997.
4. Ewald P.P. Die Berechnung optischer und elektrostatischer Gitterpotentiale // Ann. Phys. - 1921. - Vol. 64. - Р. 253-287.
5. Pereira J.C.G., Catlow C.R.A., Price G.D.J. Molecular Dynamics Simulation of Liquid H2O, MeOH, EtOH, Si(OMe)4, and Si(OEt)4, as a Function of Temperature and Pressure // Phys. Chem. - 2001. - Vol. 105. - Р. 1909-1925.
