Численное интегрирование уравнений движения

evgeshe4ka

Есть 3d задачка: в пространстве находятся частицы, между ними есть потенциал взаимодействия. Т.е. уравнений движения такие:
dx/dt = v
dv/dt = a = F(x)/m = -grad(U)/m
Нереальный тупняк - как это проинтегрировать численно?! Метод Эйлера не подходит - слишком большая неустойчивость - приходится брать слишком маленький шаг по времени. Как использовать того же Рунге-Кутта в голову совершенно не приходит...
Ну и в догонку вопрос как вставить в эту схемку моделирование теплового движения? У меня есть некоторые идейки, но они не выглядят убедительными...

KaterinKa

Да вроде используют метод Эйлера, но какой-то немного модифицированный.
Метод Рунге-Кутта в таких задачах оказывается неустойчивым.
Погугли по словам "метод молекулярной динамики".
Моделирование теплового движения производится приданием частицам в начальный момент случайных скоростей по закону mv^2/2 = 3kT/2, который должен выполняться в среднем по ансамблю.

BSCurt

 [math]$\frac{x(t_{n+1})-2x(t_n)+x(t_{n-1})}{\Delta t^2}=\frac{grad U(x(t_n}{m} $[/math] первое что прихоидт в голову.

urka3000

Чисто неявная схема для методе Эйлера абсолютно устойчива.
Почитать можно, например, у Калиткина в "Численных методах"

BSCurt

Э-э? Решается же система ОДУ, откуда ещё и шаг по х взялся.

urka3000

да, картинка не та
В общем, тут вполне уместен метод Рунге-Кутта, который называется "предиктор-корректор"

Lene81

Курить методы Адамса

evgeshe4ka

Хуже всего, что в Рунге-Кутте и Адамсе приходится считать такие вещи [math]grad $U(x_{n+1})$[/math]
[math]$x_* = x_{n+1} = x_n+v_n \Delta t+a_n \frac{\Delta^2}{2}$[/math] можно посчитать методом Эйлера в качестве промежуточного. А потом использовать его для вычисления новых ускорений и уточненного значения [math]$x_{n+1}$[/math]. Попробовал так, но результат совсем не впечатляет

sashok01

Хуже всего, что в Рунге-Кутте и Адамсе приходится считать такие вещи как [math]$\mathrm{grad}\, U(x_{n+1})$[/math]
так в этов и вся суть неявных схем. Они сложнее реализуются, зато более устойчивы.

evgeshe4ka

так в этов и вся суть неявных схем. Они сложнее реализуются, зато более устойчивы.
спасибо, кэп =) ток не обижайся )
Просто недавно одногруппник зашел, говорит: ты глупый что ли? не слушал что ли про неявные схемы?! я тебе сейчас википедию покажу - я его в дверь выпихал

demiurg

Это же молекулярная динамика. Используют чаще всего алгоритм Верле (Verlet). Это, собственно, прямо в лоб, от эйлера отличается только тем, что следующая точка по времени вычисляется не по одной, а по двум предыдущим.
http://en.wikipedia.org/wiki/Verlet_integration
А вообще лучше напиши, что тебе реально надо, может там ещё другие проблемы есть, или вообще лучше по-другому зайти.

evgeshe4ka

Я подозреваю, что что-то не в порядке с числами в потенциалах, которые я использую...
В общем, тот же Верлет использует [math]$\Delta е$[/math] порядка 1e-15 - фемтосекунды - примерно одна десятая самых быстрых молекулярных колебаний. У меня с такими временами получается, что смещения порядка ангстремов - это слишком много.
Также непонятно, как вычислять градиент. Самое простое это так [math]$ma_x = -{grad U}_x=- \frac{U(x+\Delta x)-U(x-\Delta x)}{2\Delta x}$[/math]
 [math]$\Delta x$[/math] я беру около 0.001 ангстрема - в принципе должно хватать .
Чтобы решение было хоть как-то устойчивым приходится брать [math]$\Delta t = 10^{-17}$[/math] секунды - это слишком маленький шаг. Время моделирования неприемлемо большое

evgeshe4ka

А вообще задача простая:
Несколько молекул воды могут образовывать стационарное состояние. Димер воды - существует и устойчив. Надо для больших структур построить такие устойчивые образования. И посмотреть как такие структуры будут вести себя, если учитывать тепловое движение атомов/молекул. Т.е. фактически найти температуру при которой они разваливаются.
координаты атомов задаются изначально
потенциалы взаимодействия известны
рассчитываются скорости атомов
и вычисляются положения атомов
смотрится как такая структуры ведет себя
и потом добавляется тепло kT

urka3000

По-моему, здесь не помешает нормировка.
Порядков на 13-15.

sashok01

ну а в чем разница? у тебя та же задача Коши, только не для одной переменной, а для вектора из 6N переменных.
[math]  \begin{equation*}  \begin{gathered}  \begin{bmatrix} \dot p\\\dot q\end{bmatrix}  = \begin{bmatrix} F_1(q,p)\\F_2(q,p)\end{bmatrix} = \begin{bmatrix} q\\-\mathrm{grad} \, U(p)\end{bmatrix}\\  p=\begin{bmatrix}x^1_1\\x^1_2\\x^1_3 \\ \vdots \\ x^N_1\\x^N_2\\x^N_3\end{bmatrix}, q = \dot p = \begin{bmatrix} \\ \dot x^1_1\\\dot x^1_2\\\dot x^1_3 \\ \vdots \\ \dot x^N_1\\ \dot x^N_2\\ \dot x^N_3\end{bmatrix}  \end{gathered}  \end{equation*}  [/math]
Например, неявная схема эйлера:
[math]  \begin{equation*}  \begin{gathered}  \begin{bmatrix} p_{n+1}\\ q_{n+1}\end{bmatrix}  = \begin{bmatrix}p_{n}\\ q_{n}\end{bmatrix} + \tau \begin{bmatrix} q_{n+1}\\-\mathrm{grad} \, U(p_{n+1})\end{bmatrix}\\  \end{gathered}  \end{equation*}  [/math]
Чтобы не решать систему нелинейных уравнений, вектор [math]$ \begin{bmatrix} q_{n+1}\\-\mathrm{grad} \, U(p_{n+1})\end{bmatrix}\\$[/math] находят приближенно явной схемой эйлера
PS. А почему градиент считается приближенно конечными разностями? можно же точную формулу написать, по порядку сложности разве будет намного сложнее вычисления самого потенциала? Когда я решал численно системы нелинейных уравнений, производные,взятые численно, мешались. После того, как я их заменил на точные, все стало намного лучше.

demiurg

01 фемтосекунды много, да. Должно 1-2 хватать.
А какие потенциалы? Про это же, в общем, много чего сделано и понаписано. Есть всякие там модели воды:TIP3P, SPC, TIP4P...
Ну и в догонку вопрос как вставить в эту схемку моделирование теплового движения? У меня есть некоторые идейки, но они не выглядят убедительными...
Есть несколько варивантов: термостат Nose-Hoover, термостат Berendsen или Andersen. Не очень помню чем они между собой отличаются, но непринципиально. В каком-то из них есть фиктивная частица, с которой взаимодействют все остальные, и её скорость каждый шаг меняется так, чтобы температура была постоянной. А какой-то другой рескейлит все скорости с той же целью.
Алтернативно, можно делать Ланжевеновскую (броуновскую) динамику, где на каждую частицу действует вдобавок ещё случайная сила, в зависимости от температуры. Ну это тебе не очень подходит, потому что подразумевается, что это действие растворителя, и там есть ещё параметр вязкости, который устанавливает timescale фактически.
Я бы для твоей цели вообще взял готовый пакет, например LAMMPS, там можно всё это делать. Ну или там CHARMM, AMBER или GROMACS, но они не такие гибкие.

evgeshe4ka

Чтобы не решать систему нелинейных уравнений, вектор [math]$$\begin{bmatrix} q_{n+1} \\ -gradU(p_{n+1}) \end{bmatrix} \qquad$$[/math] находят приближенно явной схемой эйлера
Спасибо огромное! А можно немного поподробнее про это. Ты хочешь сказать, что вычисляется [math]$$\begin{bmatrix} q_{n+1} \\ p_{n+1} \end{bmatrix} \qquad$$[/math] явно, а потом что?

sashok01

ну да, сначала находят [math] $\begin{bmatrix}\hat p_{n+1} \\ \hat q_{n+1}\end{bmatrix}$ [/math] приближенно явной схемой эйлера:
[math]  \begin{equation*}  \begin{bmatrix} \hat p_{n+1}\\ \hat q_{n+1}\end{bmatrix}  = \begin{bmatrix}p_{n}\\ q_{n}\end{bmatrix} + \tau \begin{bmatrix} q_{n}\\-\mathrm{grad} \, U(p_{n})\end{bmatrix}\\  \end{equation*}  [/math],
а потом подставляют в неявную схему эйлера
[math]  \begin{equation*}  \begin{bmatrix}  p_{n+1}\\ q_{n+1}\end{bmatrix}  = \begin{bmatrix}p_{n}\\ q_{n}\end{bmatrix} + \tau \begin{bmatrix} \hat q_{n+1}\\-\mathrm{grad} \, U(\hat p_{n+1})\end{bmatrix}\\  \end{equation*}  [/math]
Аналогично поступают с методом Рунге-Кутта. Я думаю, что тут это оправдано, поскольку точная система является нелинейной ( потенциал U - нелинейный скорее всего). Например, в линейной теории упругости так не делают - честно применяют конечно-разностный аналог линейного дифференциального оператора на неизвестном временном слое.

evgeshe4ka

А какие потенциалы?
Полтева-Маленкова. Смысл в том, что HyperChem и другие стандартные программы считают в закрытом виде. Используют в большинстве своем Amber. Тепло туда впихнуть очень сложно. А надо посмотреть как разные потенциалы ведут себя при учете тепла и т.д.

evgeshe4ka

Аналогично поступают с методом Рунге-Кутта
Ну я и поступил аналогично - ничего хорошего. Чует сердце, что надо всю нелинейную схему в неявном виде расписывать...

demiurg

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

demiurg

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

sslk

в namd удобно регулировать температуру и минимизировать энергию
а почему такой потенциал выбран?

sashok01

Кстати, а какие проблемы с Рунге -Кутта? Там вроде бы все естественным образом обобщается на векторное ОДУ.
[math]  \begin{equation*}  \begin{gathered}  \vec {\dot x} =\vec f(\vec x)\\  \vec x_{n+1} = \vec x_n  + \tau /6 (\vec k_1 + 2 \vec k_2 +2 \vec k_3 + \vec k_4) \\  \vec k_1=\vec f(\vec x_n)\\  \vec k_2=\vec f(\vec x_n+\tau/2\vec k_1)\\  \vec k_3=\vec f(\vec x_n+\tau/2\vec k_2)\\  \vec k_4=\vec f(\vec x_n+\tau \vec k_3)  \end{gathered}  \end{equation*}  [/math]

Brina

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

Brina

Полтева-Маленкова
Можешь формулу написать?

KaterinKa

Да вроде говорят, что метод Рунге-Кутта неустойчив на больших временах и для большого числа частиц.
Чем больше точность, тем меньше устойчивость.

Brina

Тогда, наверное, как уже предлагали использовать неявную схему. Она железобетонная, проблема может вылезти в численном решении алгебраических уравнений... А как пофиксить эту тему — ХЗ...
Есть еще одна тривиальная идейка. Если прога сильно завязана на скорость, то ее можно распараллелить — такие задачи, насколько я себе представляю, неплохо параллелятся при разбиения области ингегрирования по пространству (хотя проблема может вылезти, когда большая часть частиц соберется в одной области — тогда распараллеливание теряет свою ценность).

KaterinKa

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

Brina

Это да... Согласен абсолютно.

Brina

И заголовок "Численное интегрирование уравнений движения" не совсем соотвестствует теме...

terl

А вообще задача простая:
я думаю уместно было бы начать работу над этой проблемой с создания ЭВМ!

BSCurt

я думаю уместно было бы начать работу над этой проблемой с создания ЭВМ!
Не читери начинать надо с колеса и далее по TechTree.

andriano_as

Мы считали недавно теплопроводность нанотрубок и графена с помощью нашей молдинамической программы MD-kMC методом неравновесной молекулярной динамики (обратный метод NEMD термостаты Берендсена использовались для поддержания теплового потока. Общая схема расчета выглядела так:
.
Интегрирование уравнений движения в MD-kMC выполняется, кажется, методом Верле. Одной из основных проблем у нас было достаточно большое время установления стационарного теплового потока в длинных листах графена. Расчеты для узкого, но длинного листа графена (1 мкм) выполнялись не меньше недели на 200 ядрах. Для моделирования межатомного взаимодействия использовался упрощенный потенциал Бреннера. Кое-что из результатов есть здесь.

andriano_as

Да, еще есть метод Бимана, если хочется самому возиться с численными схемами. Вообще по методам и программам лучше почитать эту статью.
Оставить комментарий
Имя или ник:
Комментарий: