Адаптация алгоритма Верле под многопроцессорные системы Текст научной статьи по специальности «Физика»

2. Информационные системы и технологии

УДК 519.688

© Е.И. Герман, Ш.Б. Цыдыпов

АДАПТАЦИЯ АЛГОРИТМА ВЕРЛЕ ПОД МНОГОПРОЦЕССОРНЫЕ СИСТЕМЫ

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

Ключевые слова: параллельные вычисления, метод молекулярной динамики, алгоритм Верле.

E.I. Herman, Sh.B. Tsydypov

ADAPTATION OF ALGORITHM VERLET FOR MULTIPROCESSOR SYSTEMS

The methodology of implementation of parallel asynchronous computations is expounded for modeling molecular systems using molecular dynamics simulations.

Keywords: parallel computing, molecular dynamics method, the Verlet algorithm.

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

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

Метод молекулярной динамики, попросту говоря, является численной реализацией решения уравнений движения Ньютона для множества частиц. Полную силу, действующую на i-й атом, можно представить в форме суммы векторов [2]:

N-1 , ч

F =-Xv°(i/). (1)

При этом предполагается, что известны координаты центров всех атомов и вид потенциала взаимодействия Ф(г). Энергия частиц инертного газа, например, может определяться потенциалом Леннарда-Джонса:

ф(г) = 4s

i r -, 12 r -, 6 Л

a a

r r

V /

(2)

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

1 N-1

д1 (? + Д) = ц «)Д - —£УФ(у)(Д?)2 + ($), (3)

2т,, у

1 N-1

+ Д) = ц(0--£УФ(//)Д . (4)

i J

Возьмем для примера систему из N частиц инертного газа.

Алгоритм Верле позволяет повысить точность итерации моделирования за счет частичного изменения скоростей частиц [3]. В общем случае одна итерация данного алгоритма на языке C++ представляет вид: void Verlet(void){

for (int i=0; i!=N; i++){

/* Расчет координаты i-й частицы*/ particles [i].x = particles[i].x + particles[i].vx*dt +

particle s [i]. ax* dt2_2; particles[i].y = particles[i].y + particles[i].vy*dt +

particles[i].ay* dt2_2; particles[i].z = particles [i].z + particles[i].vz*dt +

particle s[i].az* dt2_2; /* Корректировка координат частицы согласно условий

периодических граничных условий*/ periodic(particles[i].x, particles[i].y, particles[i].z); /*Частичное изменение скорости */ particles [i]. vx+=particles[i].ax*dt*0.5; particle s [i]. vy+=particles [i] .ay * dt* 0. 5; particles [i] .vz+=particles [i] .az* dt* 0.5;

}

GetAccelO^/расчет ускорений for (int i=0; i!=N; i++){

/*Изменение скорости с учетом вычисленного ускорения*/ particles [i]. vx+=particles[i].ax*dt*0.5; particle s [i]. vy+=particles [i]. ay * dt* 0. 5; particles [i] .vz+=particles [i] .az* dt* 0.5;

}

}

В представленном листинге производится вызов подпрограммы GetAccel(), задачей которой является определение ускорений частиц ис-

ходя из (4). Видно, что для определения сил, действующих на одну частицу, необходимо определять расстояния до всех остальных N-1 частиц. Эта процедура дает наибольшую вычислительную нагрузку. В классическом представлении на языке C++ код подпрограммы выглядит следующим образом:

void GetAccel(void){

for (int i=0; i!=N; i++){ particles[i].ax=ü; particles[i].ay=0; particles[i].az=ü;

}

for (int i=0; i!=N-1; i++){

for (int j=i+1; j!=N; j++){

dx=particles [j ] .x-particles [i] .x; dy=particles[j] .y-particles[i] .y; dz=particles [j ]. z-particle s [i]. z; Separate(dx,dy,dz); r=sqrt(dx* dx+dy * dy+dz* dz); r=1/r; r2=r*r; r3=r2*r; r6=r3*r3; a=24*r2*r6*(2*r6-1); particles [i] .ax+=a*dx; particles [i] .ay+=a* dy; particles [i] .az+=a* dz; particles [j ] .ax-=a* dx; particles [j ] .ay-=a* dy; particles [j ] .az-=a* dz;

}

}

}

В среде программирования MS Visual Studio имеется компонент Net. FrameworkBackgroundWorker, с помощью которого возможна организация вычислительных операций в асинхронном фоновом режиме.

Массив координат частиц системы при расчете ускорений представляет собой массив связанных переменных, здесь ускорение i-й частицы напрямую связано с положением любой другой j-й частицы, поэтому прямое разделение подпрограммы Get Accel () на два и более асинхронных процесса фактически невозможно. Вместе с тем возможно использование одного управляющего компонента BackgroundWorker0 и двух счетных компонентов BackgroundWorker1 и BackgroundWorker2. Тогда схема разделения вычислительных операций одной итерации МД будет представлять следующую последовательность:

1. BackgroundWorkerü обнуляет ускорения на начале данной итерации, запускает BackgroundWorker1 и BackgroundWorker2 и ожидает событий завершения работы этих компонентов.

2. С использованием компонента BackgroundWorker1 производится расчет ускорений по первым N/2 частицам системы, с использованием BackgroundWorker2 - по последним N/2 частицам.

3. По завершении выполнения работы асинхронных потоков при помощи BackgroundWorker0 производится расчет ускорений, вызванных силами взаимодействия первых N/2 частиц и последних N/2 частиц.

Рис 1. Схема разделения подпрограммы ве1Лссе1() на вычислительные потоки

Представленная схема использована для двухпроцессорной вычислительной системы. Теоретически прямой расчет ускорений для системы из тысячи частиц влечет за собой 499500 итераций цикла, при использовании предложенной схемы разделения вычислений на один процессор приходится 249500 итераций, что приводит к сокращению времени расчетов примерно в два раза. На практике увеличение скорости вычислений составляет порядка 60%. Это связано со спецификой работы операционной системы и неравнозначностью вычислительных ядер процессора.

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

Литература

1. Герман Е.И., Цыдыпов Ш.Б. Радиальные функции распределения неравновесных систем // Вестник Бурятского государственного университета. - 2013. - Вып. 3. - С. 104-107.

2. Хеерман Д.В. Методы компьютерного эксперимента в физике. - М.: Наука, 1990.

3. Гулд. Х., Тобочник Я. Компьютерное моделирование в физике. - М.: Мир, 1990.

Цыдыпов Шулун Балдоржиевич, доктор технических наук, доцент, заведующий кафедрой общей физики физико-технического факультета Бурятского госуниверситета. Тел. 89025634889, e-mail: shulun@bsu.ru

Герман Евгений Иванович, преподаватель кафедры общей физики физико-технического факультета Бурятского госуниверситета. Тел. 89244560177, e-mail: net-admin@list.ru

Tsydypov Shulun Baldorzhievich, doctor of technical sciences, associate professor, head of the department of general physics, Buryat State University. Tel: 89025634889, e-mail: shulun@bsu.ru

Herman Evgeny Ivanovich, lecturer, department of general physics, Buryat State University. Tel: 89244560177, e-mail: net-admin@list.ru