Алгоритмы интегрирования уравнений движения в приложении метода молекулярной динамики

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

ММД основывается на численном решении систем дифференциальных уравнений движения вида:

d2

dt2

d2

dt2

mt

i(

=1..N у

j*t

i)

—»

Ш;

tj

f t

=1..N у

j*t

Для численного решения необходимо представить уравнение (2), с учетом равенства (3), в виде системы уравнений первого порядка.

d^,

—>

-77 = f dt

d"rt

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

Для итерационной схемы верны следующие равенства:

't^ n+1 = '^^(t + dt)                                                           (5)

i^n+i = i^(t + dt)                                                             ⁽⁶⁾

Далее, разлагая (5) и (6) в ряд Тейлора, можно получить выражения для всех нижеприведенных схем.

Начнем рассмотрение с самой простой схемы интегрирования, а именно с алгоритма Эйлера.

'^,_n+1 = ^ n + cindLt + O(dt)                                            (7)

i n+i = i n + ^ n dt + O(dt)                                        (8)

Как показано в (7,8), алгоритм Эйлера имеет первый порядок точности по координатам и скоростям. Главным недостатком данной схемы является асимметричность относительно шага интегрирования, т.е. при продвижении решения на временной шаг dt используется информация о производных только в начальной точке интервала. Данная схема приведена, преимущественно, в качестве примера, так как на практике не используется.

Следующим рассмотрим алгоритм Эйлера-Кромера.

^ n+i = ^ n + d ^ dt + O(dt)

i n+i = ’ + ^ n+i dt + O(dt)

Хотя схема (9,10) имеет, аналогично (7,8), первый порядок точности, она является симметричной, что способствует меньшему накоплению погрешности. Данный алгоритм может быть использован в демонстрационных целях, например в рамках лабораторных или курсовых студенческих работ.

Последней схемой первого порядка точности является метод центральной точки.

4" _п+1 = 4"П + didt + O(dt)                                         (11)

i n+i = i n + № + ^ n+i)^ + O(dt 2 )                           ⁽¹²⁾

В отличие от алгоритмов Эйлера и Эйлера-Кромера, данный метод имеет второй порядок точности по координате.

Перейдем к рассмотрению более точных схем интегрирования. Наиболее простой и распространенной из них является алгоритм Верле в скоростной форме.

i_n+i = i n + ^ n dt + -C_ndt² + O(dt³)

^" n+i    ^"n ⁺ 2 ^(c n+i ⁺ dn)^{dt + O(dt} )

Основным достоинством данной схемы является простота реализации. Как можно заметить из (14), для корректной работы необходимо сохранять в памяти дополнительные значения ускорений. Данная схема может применяться для реализации численных экспериментов, направленных на выявление статистических или макроскопических явлений.

Последней рассмотренной схемой является алгоритм Бимана в явной форме.

i^t_n+i = i n + ^;i’_ndt + — (4ct_n — 6i_n-i)dt² + O(dt⁴)                      (15)

^ _п+1 = ^ _п + 7 (2c_n+i + 5c"_n — c"_n-i)dt + O(dt³)                   ⁽¹⁶⁾

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