Численное интегрирование уравнений движения
Да вроде используют метод Эйлера, но какой-то немного модифицированный.
Метод Рунге-Кутта в таких задачах оказывается неустойчивым.
Погугли по словам "метод молекулярной динамики".
Моделирование теплового движения производится приданием частицам в начальный момент случайных скоростей по закону mv^2/2 = 3kT/2, который должен выполняться в среднем по ансамблю.
Метод Рунге-Кутта в таких задачах оказывается неустойчивым.
Погугли по словам "метод молекулярной динамики".
Моделирование теплового движения производится приданием частицам в начальный момент случайных скоростей по закону mv^2/2 = 3kT/2, который должен выполняться в среднем по ансамблю.
Чисто неявная схема для методе Эйлера абсолютно устойчива.
Почитать можно, например, у Калиткина в "Численных методах"
Почитать можно, например, у Калиткина в "Численных методах"
Э-э? Решается же система ОДУ, откуда ещё и шаг по х взялся.
да, картинка не та
В общем, тут вполне уместен метод Рунге-Кутта, который называется "предиктор-корректор"
В общем, тут вполне уместен метод Рунге-Кутта, который называется "предиктор-корректор"
Курить методы Адамса
Хуже всего, что в Рунге-Кутте и Адамсе приходится считать такие вещи ![[math]grad $U(x_{n+1})$[/math]](mathimg.php?math=grad%20%24U%28x_%7Bn%2B1%7D%29%24)
![[math]$x_* = x_{n+1} = x_n+v_n \Delta t+a_n \frac{\Delta^2}{2}$[/math]](mathimg.php?math=%24x_%2A%20%3D%20x_%7Bn%2B1%7D%20%3D%20x_n%2Bv_n%20%5CDelta%20t%2Ba_n%20%5Cfrac%7B%5CDelta%5E2%7D%7B2%7D%24)
можно посчитать методом Эйлера в качестве промежуточного. А потом использовать его для вычисления новых ускорений и уточненного значения ![[math]$x_{n+1}$[/math]](mathimg.php?math=%24x_%7Bn%2B1%7D%24)
. Попробовал так, но результат совсем не впечатляет
Хуже всего, что в Рунге-Кутте и Адамсе приходится считать такие вещи кактак в этов и вся суть неявных схем. Они сложнее реализуются, зато более устойчивы.
так в этов и вся суть неявных схем. Они сложнее реализуются, зато более устойчивы.спасибо, кэп =) ток не обижайся )
Просто недавно одногруппник зашел, говорит: ты глупый что ли? не слушал что ли про неявные схемы?! я тебе сейчас википедию покажу - я его в дверь выпихал
Это же молекулярная динамика. Используют чаще всего алгоритм Верле (Verlet). Это, собственно, прямо в лоб, от эйлера отличается только тем, что следующая точка по времени вычисляется не по одной, а по двум предыдущим.
http://en.wikipedia.org/wiki/Verlet_integration
А вообще лучше напиши, что тебе реально надо, может там ещё другие проблемы есть, или вообще лучше по-другому зайти.
http://en.wikipedia.org/wiki/Verlet_integration
А вообще лучше напиши, что тебе реально надо, может там ещё другие проблемы есть, или вообще лучше по-другому зайти.
Я подозреваю, что что-то не в порядке с числами в потенциалах, которые я использую...
В общем, тот же Верлет использует![[math]$\Delta е$[/math]](mathimg.php?math=%24%5CDelta%20%E5%24)
порядка 1e-15 - фемтосекунды - примерно одна десятая самых быстрых молекулярных колебаний. У меня с такими временами получается, что смещения порядка ангстремов - это слишком много.
Также непонятно, как вычислять градиент. Самое простое это так![[math]$ma_x = -{grad U}_x=- \frac{U(x+\Delta x)-U(x-\Delta x)}{2\Delta x}$[/math]](mathimg.php?math=%24ma_x%20%3D%20-%7Bgrad%20U%7D_x%3D-%20%5Cfrac%7BU%28x%2B%5CDelta%20x%29-U%28x-%5CDelta%20x%29%7D%7B2%5CDelta%20x%7D%24)
![[math]$\Delta x$[/math]](mathimg.php?math=%24%5CDelta%20x%24)
я беру около 0.001 ангстрема - в принципе должно хватать .
Чтобы решение было хоть как-то устойчивым приходится брать![[math]$\Delta t = 10^{-17}$[/math]](mathimg.php?math=%24%5CDelta%20t%20%3D%2010%5E%7B-17%7D%24)
секунды - это слишком маленький шаг. Время моделирования неприемлемо большое
В общем, тот же Верлет использует
Также непонятно, как вычислять градиент. Самое простое это так
Чтобы решение было хоть как-то устойчивым приходится брать
А вообще задача простая:
Несколько молекул воды могут образовывать стационарное состояние. Димер воды - существует и устойчив. Надо для больших структур построить такие устойчивые образования. И посмотреть как такие структуры будут вести себя, если учитывать тепловое движение атомов/молекул. Т.е. фактически найти температуру при которой они разваливаются.
координаты атомов задаются изначально
потенциалы взаимодействия известны
рассчитываются скорости атомов
и вычисляются положения атомов
смотрится как такая структуры ведет себя
и потом добавляется тепло kT
Несколько молекул воды могут образовывать стационарное состояние. Димер воды - существует и устойчив. Надо для больших структур построить такие устойчивые образования. И посмотреть как такие структуры будут вести себя, если учитывать тепловое движение атомов/молекул. Т.е. фактически найти температуру при которой они разваливаются.
координаты атомов задаются изначально
потенциалы взаимодействия известны
рассчитываются скорости атомов
и вычисляются положения атомов
смотрится как такая структуры ведет себя
и потом добавляется тепло kT
По-моему, здесь не помешает нормировка.
Порядков на 13-15.
Порядков на 13-15.
ну а в чем разница? у тебя та же задача Коши, только не для одной переменной, а для вектора из 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]](mathimg.php?math=%0D%0A%5Cbegin%7Bequation%2A%7D%0D%0A%5Cbegin%7Bgathered%7D%0D%0A%5Cbegin%7Bbmatrix%7D%20%5Cdot%20p%5C%5C%5Cdot%20q%5Cend%7Bbmatrix%7D%20%20%3D%20%5Cbegin%7Bbmatrix%7D%20F_1%28q%2Cp%29%5C%5CF_2%28q%2Cp%29%5Cend%7Bbmatrix%7D%20%3D%20%5Cbegin%7Bbmatrix%7D%20q%5C%5C-%5Cmathrm%7Bgrad%7D%20%5C%2C%20U%28p%29%5Cend%7Bbmatrix%7D%5C%5C%0D%0Ap%3D%5Cbegin%7Bbmatrix%7Dx%5E1_1%5C%5Cx%5E1_2%5C%5Cx%5E1_3%20%5C%5C%20%5Cvdots%20%5C%5C%20x%5EN_1%5C%5Cx%5EN_2%5C%5Cx%5EN_3%5Cend%7Bbmatrix%7D%2C%20q%20%3D%20%5Cdot%20p%20%3D%20%5Cbegin%7Bbmatrix%7D%20%5C%5C%20%5Cdot%20x%5E1_1%5C%5C%5Cdot%20x%5E1_2%5C%5C%5Cdot%20x%5E1_3%20%5C%5C%20%5Cvdots%20%5C%5C%20%5Cdot%20x%5EN_1%5C%5C%20%5Cdot%20x%5EN_2%5C%5C%20%5Cdot%20x%5EN_3%5Cend%7Bbmatrix%7D%0D%0A%5Cend%7Bgathered%7D%0D%0A%5Cend%7Bequation%2A%7D%0D%0A)
Например, неявная схема эйлера:
![[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]](mathimg.php?math=%0D%0A%5Cbegin%7Bequation%2A%7D%0D%0A%5Cbegin%7Bgathered%7D%0D%0A%5Cbegin%7Bbmatrix%7D%20p_%7Bn%2B1%7D%5C%5C%20q_%7Bn%2B1%7D%5Cend%7Bbmatrix%7D%20%20%3D%20%5Cbegin%7Bbmatrix%7Dp_%7Bn%7D%5C%5C%20q_%7Bn%7D%5Cend%7Bbmatrix%7D%20%2B%20%5Ctau%20%5Cbegin%7Bbmatrix%7D%20q_%7Bn%2B1%7D%5C%5C-%5Cmathrm%7Bgrad%7D%20%5C%2C%20U%28p_%7Bn%2B1%7D%29%5Cend%7Bbmatrix%7D%5C%5C%0D%0A%5Cend%7Bgathered%7D%0D%0A%5Cend%7Bequation%2A%7D%0D%0A)
Чтобы не решать систему нелинейных уравнений, вектор![[math]$ \begin{bmatrix} q_{n+1}\\-\mathrm{grad} \, U(p_{n+1})\end{bmatrix}\\$[/math]](mathimg.php?math=%24%20%5Cbegin%7Bbmatrix%7D%20q_%7Bn%2B1%7D%5C%5C-%5Cmathrm%7Bgrad%7D%20%5C%2C%20U%28p_%7Bn%2B1%7D%29%5Cend%7Bbmatrix%7D%5C%5C%24)
находят приближенно явной схемой эйлера
PS. А почему градиент считается приближенно конечными разностями? можно же точную формулу написать, по порядку сложности разве будет намного сложнее вычисления самого потенциала? Когда я решал численно системы нелинейных уравнений, производные,взятые численно, мешались. После того, как я их заменил на точные, все стало намного лучше.
Например, неявная схема эйлера:
Чтобы не решать систему нелинейных уравнений, вектор
PS. А почему градиент считается приближенно конечными разностями? можно же точную формулу написать, по порядку сложности разве будет намного сложнее вычисления самого потенциала? Когда я решал численно системы нелинейных уравнений, производные,взятые численно, мешались. После того, как я их заменил на точные, все стало намного лучше.
01 фемтосекунды много, да. Должно 1-2 хватать.
А какие потенциалы? Про это же, в общем, много чего сделано и понаписано. Есть всякие там модели воды:TIP3P, SPC, TIP4P...
Алтернативно, можно делать Ланжевеновскую (броуновскую) динамику, где на каждую частицу действует вдобавок ещё случайная сила, в зависимости от температуры. Ну это тебе не очень подходит, потому что подразумевается, что это действие растворителя, и там есть ещё параметр вязкости, который устанавливает timescale фактически.
Я бы для твоей цели вообще взял готовый пакет, например LAMMPS, там можно всё это делать. Ну или там CHARMM, AMBER или GROMACS, но они не такие гибкие.
А какие потенциалы? Про это же, в общем, много чего сделано и понаписано. Есть всякие там модели воды:TIP3P, SPC, TIP4P...
Ну и в догонку вопрос как вставить в эту схемку моделирование теплового движения? У меня есть некоторые идейки, но они не выглядят убедительными...Есть несколько варивантов: термостат Nose-Hoover, термостат Berendsen или Andersen. Не очень помню чем они между собой отличаются, но непринципиально. В каком-то из них есть фиктивная частица, с которой взаимодействют все остальные, и её скорость каждый шаг меняется так, чтобы температура была постоянной. А какой-то другой рескейлит все скорости с той же целью.
Алтернативно, можно делать Ланжевеновскую (броуновскую) динамику, где на каждую частицу действует вдобавок ещё случайная сила, в зависимости от температуры. Ну это тебе не очень подходит, потому что подразумевается, что это действие растворителя, и там есть ещё параметр вязкости, который устанавливает timescale фактически.
Я бы для твоей цели вообще взял готовый пакет, например LAMMPS, там можно всё это делать. Ну или там CHARMM, AMBER или GROMACS, но они не такие гибкие.
Чтобы не решать систему нелинейных уравнений, векторСпасибо огромное! А можно немного поподробнее про это. Ты хочешь сказать, что вычисляетсянаходят приближенно явной схемой эйлера
ну да, сначала находят ![[math] $\begin{bmatrix}\hat p_{n+1} \\ \hat q_{n+1}\end{bmatrix}$ [/math]](mathimg.php?math=%20%24%5Cbegin%7Bbmatrix%7D%5Chat%20p_%7Bn%2B1%7D%20%5C%5C%20%5Chat%20q_%7Bn%2B1%7D%5Cend%7Bbmatrix%7D%24%20)
приближенно явной схемой эйлера:
![[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]](mathimg.php?math=%0D%0A%5Cbegin%7Bequation%2A%7D%0D%0A%5Cbegin%7Bbmatrix%7D%20%5Chat%20p_%7Bn%2B1%7D%5C%5C%20%5Chat%20q_%7Bn%2B1%7D%5Cend%7Bbmatrix%7D%20%20%3D%20%5Cbegin%7Bbmatrix%7Dp_%7Bn%7D%5C%5C%20q_%7Bn%7D%5Cend%7Bbmatrix%7D%20%2B%20%5Ctau%20%5Cbegin%7Bbmatrix%7D%20q_%7Bn%7D%5C%5C-%5Cmathrm%7Bgrad%7D%20%5C%2C%20U%28p_%7Bn%7D%29%5Cend%7Bbmatrix%7D%5C%5C%0D%0A%5Cend%7Bequation%2A%7D%0D%0A)
,
а потом подставляют в неявную схему эйлера
![[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]](mathimg.php?math=%0D%0A%5Cbegin%7Bequation%2A%7D%0D%0A%5Cbegin%7Bbmatrix%7D%20%20p_%7Bn%2B1%7D%5C%5C%20q_%7Bn%2B1%7D%5Cend%7Bbmatrix%7D%20%20%3D%20%5Cbegin%7Bbmatrix%7Dp_%7Bn%7D%5C%5C%20q_%7Bn%7D%5Cend%7Bbmatrix%7D%20%2B%20%5Ctau%20%5Cbegin%7Bbmatrix%7D%20%5Chat%20q_%7Bn%2B1%7D%5C%5C-%5Cmathrm%7Bgrad%7D%20%5C%2C%20U%28%5Chat%20p_%7Bn%2B1%7D%29%5Cend%7Bbmatrix%7D%5C%5C%0D%0A%5Cend%7Bequation%2A%7D%0D%0A)
Аналогично поступают с методом Рунге-Кутта. Я думаю, что тут это оправдано, поскольку точная система является нелинейной ( потенциал U - нелинейный скорее всего). Например, в линейной теории упругости так не делают - честно применяют конечно-разностный аналог линейного дифференциального оператора на неизвестном временном слое.
а потом подставляют в неявную схему эйлера
Аналогично поступают с методом Рунге-Кутта. Я думаю, что тут это оправдано, поскольку точная система является нелинейной ( потенциал U - нелинейный скорее всего). Например, в линейной теории упругости так не делают - честно применяют конечно-разностный аналог линейного дифференциального оператора на неизвестном временном слое.
А какие потенциалы?Полтева-Маленкова. Смысл в том, что HyperChem и другие стандартные программы считают в закрытом виде. Используют в большинстве своем Amber. Тепло туда впихнуть очень сложно. А надо посмотреть как разные потенциалы ведут себя при учете тепла и т.д.
Аналогично поступают с методом Рунге-КуттаНу я и поступил аналогично - ничего хорошего. Чует сердце, что надо всю нелинейную схему в неявном виде расписывать...
LAMMPS в открытом. В амбере, кстати, термостаты тоже есть.
Но я к тому, что у тебя дело в ошибке, скорее всего, а не в том, что плохая разностная схема. Схема верле у всех сходится. Не должно быть смещений на ангстремы за фемтосекунду, ну никак. Это на первом шаге так у тебя?
Но я к тому, что у тебя дело в ошибке, скорее всего, а не в том, что плохая разностная схема. Схема верле у всех сходится. Не должно быть смещений на ангстремы за фемтосекунду, ну никак. Это на первом шаге так у тебя?
Может ещё ты их совсем неправильно ставишь вначале, молекулы в смысле. Поэтому большие силы получаются. Прежде чем гонять динамику можно сделать просто минимизацию (в смысле без времени и скоростей).
в namd удобно регулировать температуру и минимизировать энергию
а почему такой потенциал выбран?
а почему такой потенциал выбран?
Кстати, а какие проблемы с Рунге -Кутта? Там вроде бы все естественным образом обобщается на векторное ОДУ.
![[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]](mathimg.php?math=%0D%0A%5Cbegin%7Bequation%2A%7D%0D%0A%5Cbegin%7Bgathered%7D%0D%0A%5Cvec%20%7B%5Cdot%20x%7D%20%3D%5Cvec%20f%28%5Cvec%20x%29%5C%5C%0D%0A%5Cvec%20x_%7Bn%2B1%7D%20%3D%20%5Cvec%20x_n%20%20%2B%20%5Ctau%20%2F6%20%28%5Cvec%20k_1%20%2B%202%20%5Cvec%20k_2%20%2B2%20%5Cvec%20k_3%20%2B%20%5Cvec%20k_4%29%20%5C%5C%0D%0A%5Cvec%20k_1%3D%5Cvec%20f%28%5Cvec%20x_n%29%5C%5C%0D%0A%5Cvec%20k_2%3D%5Cvec%20f%28%5Cvec%20x_n%2B%5Ctau%2F2%5Cvec%20k_1%29%5C%5C%0D%0A%5Cvec%20k_3%3D%5Cvec%20f%28%5Cvec%20x_n%2B%5Ctau%2F2%5Cvec%20k_2%29%5C%5C%0D%0A%5Cvec%20k_4%3D%5Cvec%20f%28%5Cvec%20x_n%2B%5Ctau%20%5Cvec%20k_3%29%0D%0A%5Cend%7Bgathered%7D%0D%0A%5Cend%7Bequation%2A%7D%0D%0A)
Я тоже не вижу проблем с методом Р-К для систем обыкновенных уравнений. На втором курсе физфака на проге решается много задач такого рода — типа движения в каком-то потенциале.
Я бы честно попробовал выписать потенциальную энергию U, а потом бы попробовал и посчитать градиент — его расчет сводится все таки к взятию производных (причем первых чему, насколько мне известно, можно и зайца научить.
Я бы честно попробовал выписать потенциальную энергию U, а потом бы попробовал и посчитать градиент — его расчет сводится все таки к взятию производных (причем первых чему, насколько мне известно, можно и зайца научить.
Полтева-МаленковаМожешь формулу написать?
Да вроде говорят, что метод Рунге-Кутта неустойчив на больших временах и для большого числа частиц.
Чем больше точность, тем меньше устойчивость.
Чем больше точность, тем меньше устойчивость.
Тогда, наверное, как уже предлагали использовать неявную схему. Она железобетонная, проблема может вылезти в численном решении алгебраических уравнений... А как пофиксить эту тему — ХЗ...
Есть еще одна тривиальная идейка. Если прога сильно завязана на скорость, то ее можно распараллелить — такие задачи, насколько я себе представляю, неплохо параллелятся при разбиения области ингегрирования по пространству (хотя проблема может вылезти, когда большая часть частиц соберется в одной области — тогда распараллеливание теряет свою ценность).
Есть еще одна тривиальная идейка. Если прога сильно завязана на скорость, то ее можно распараллелить — такие задачи, насколько я себе представляю, неплохо параллелятся при разбиения области ингегрирования по пространству (хотя проблема может вылезти, когда большая часть частиц соберется в одной области — тогда распараллеливание теряет свою ценность).
Да молекулярная динамика — это сейчас целая наука. Для нее есть туева хуча методов, приемов и стандартных пакетов.
Пытаться использовать здесь какой-то самодельный примитивный метод — это изобретение велосипеда.
Пытаться использовать здесь какой-то самодельный примитивный метод — это изобретение велосипеда.
Это да... Согласен абсолютно.
И заголовок "Численное интегрирование уравнений движения" не совсем соотвестствует теме...
А вообще задача простая:я думаю уместно было бы начать работу над этой проблемой с создания ЭВМ!
я думаю уместно было бы начать работу над этой проблемой с создания ЭВМ!Не читери начинать надо с колеса и далее по TechTree.
Мы считали недавно теплопроводность нанотрубок и графена с помощью нашей молдинамической программы MD-kMC методом неравновесной молекулярной динамики (обратный метод NEMD термостаты Берендсена использовались для поддержания теплового потока. Общая схема расчета выглядела так:
.
Интегрирование уравнений движения в MD-kMC выполняется, кажется, методом Верле. Одной из основных проблем у нас было достаточно большое время установления стационарного теплового потока в длинных листах графена. Расчеты для узкого, но длинного листа графена (1 мкм) выполнялись не меньше недели на 200 ядрах. Для моделирования межатомного взаимодействия использовался упрощенный потенциал Бреннера. Кое-что из результатов есть здесь.
.Интегрирование уравнений движения в MD-kMC выполняется, кажется, методом Верле. Одной из основных проблем у нас было достаточно большое время установления стационарного теплового потока в длинных листах графена. Расчеты для узкого, но длинного листа графена (1 мкм) выполнялись не меньше недели на 200 ядрах. Для моделирования межатомного взаимодействия использовался упрощенный потенциал Бреннера. Кое-что из результатов есть здесь.
Да, еще есть метод Бимана, если хочется самому возиться с численными схемами. Вообще по методам и программам лучше почитать эту статью.
evgeshe4ka
Есть 3d задачка: в пространстве находятся частицы, между ними есть потенциал взаимодействия. Т.е. уравнений движения такие:dx/dt = v
dv/dt = a = F(x)/m = -grad(U)/m
Нереальный тупняк - как это проинтегрировать численно?! Метод Эйлера не подходит - слишком большая неустойчивость - приходится брать слишком маленький шаг по времени. Как использовать того же Рунге-Кутта в голову совершенно не приходит...
Ну и в догонку вопрос как вставить в эту схемку моделирование теплового движения? У меня есть некоторые идейки, но они не выглядят убедительными...