Молекулярная динамика

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

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

Варианты заданий
Расчёт производится для N молекул, где N выбирается исходя из необходимости произвести расчет для нескольких тысяч шагов по времени за приемлемый период.

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

В классической молекулярной динамике положения и скорости каждой молекулы вычисляются исходя из решения системы уравнений движения
$\frac{\large\mathbf{r}_i}{\large dt}=\mathbf{v}_i\,,\hspace{30mm}$ $\frac{\large\mathbf{v}_i}{\large dt}=\frac{\large\mathbf{F}_i}{\large m_i}\,,\hspace{30mm}$ $\mathbf{F}_i=-\frac{\large\partial U_i}{\large\partial\mathbf{r}_i}\,,$
где $i$ принимает значения от 1 до $N$ (число молекул);
$m_i\,,\hspace{10mm}\mathbf{r}_i\,,\hspace{10mm}\mathbf{v}_i$ — соответственно, масса $i$ -й частицы, положение её центра масс и скорость;
$\mathbf{F}_i$ — сумма сил, действущих на данную молекулу со стороны остальных молекул; рассчитывается исходя из суммарного потенциала $U_i=\sum\limits_{j \ne i} U_{ij}$ .

Потенциал взаимодействия двух молекул аргона удобно рассчитывать по формуле, предложенной Леннард—Джонсом:
$U_{ij}=4\varepsilon\left[\left(\frac{\large\sigma}{\large r_{ij}}\right)^{12}-\left(\frac{\large\sigma}{\large r_{ij}}\right)^6\right]\,$
Здесь $\varepsilon$ — глубина потенциальной ямы, а $\sigma$ — расстояние, на котором энергия взаимодействия равна нулю.

На практике удобно считать, что на некотором расстоянии $r_{cutoff}=\ell\varepsilon$ потенциал становится настолько мал, что при расчёте его можно строго приравнять нулю и для молекул, находящихся от данной за пределами радиуса $r_{cutoff}$ , расчёт сил взаимодействия можно не производить, что ускоряет вычисления. Величина коэффициента $\ell$ обычно выбирается из диапазона [2.5; 8].

За исключением случаев, когда возникает необходимость вычисления потенциальной или полной энергии системы, нет необходимости тратить время на непосредственное вычисление потенциала и последующее его дифференцирование. Можно произвести соответствующее алгебраическое действие и сразу считать парные силы взаимодействия
$\mathbf{f}_{ij}={\frac{\large48\varepsilon}{\large\sigma^2}}\left[\left(\frac{\large\sigma}{\large r_{ij}}\right)^{14}-\frac12\left(\frac{\large\sigma}{\large r_{ij}}\right)^8\right]\mathbf{r}_{ij}\,,$ $\hspace{30mm}\mathbf{F}_i=\sum_{j=1,\ j\ne i}^N \mathbf{f}_{ij}$

Двухэтапный метод скоростей Верле.
Для решения уравнений движения молекулярной динамики предлагается воспользоваться методом Верле в двухэтапной форме.
На первом этапе вычисляются «предварительные» знаения скоростей молекул на полушаге по времени
$\mathbf{v}_i\left(t+\frac{\large\Delta t}{\large 2}\right)=\mathbf{v}_i(t)+\frac{\large\mathbf{F}_i(t)}{\large m_i}\cdot\frac{\large\Delta t}{\large 2}$
Полученные скорости используются для вычисления новых координат молекул
$\mathbf{r}_i(t+\Delta t)=\mathbf{r}_i(t) + \mathbf{v}_i\left(t+\frac{\large\Delta t}{\large 2}\right) \cdot \Delta t$
Наконец, с учётом изменения сил, действующих на молекулы, производится окончательный пересчёт коростей
$\mathbf{v}_i(t+\Delta t)=\mathbf{v}_i\left(t+\frac{\large\Delta t}{\large 2}\right)+\frac{\large\mathbf{F}_i(t + \Delta t)}{\large m_i}\cdot\frac{\large\Delta t}{\large 2}$

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

Энергия системы.
Потенциальная:
$E_{pot}=\sum\limits_{i=1}^N U_i$
Кинетическая:
$E_{kin}=\sum\limits_{i=1}^N \frac{\large m_i |\mathbf{v}_i|^2}{\large 2}$
Полная:
$E=E_{pot} + E_{kin}$
Температура системы.
Для системы одноатомных молекул температура рассчитывается как
$T=\frac{\large 2}{\large 3k_b}E_{kin}$
где $k_b$ — константа Больцмана.
Давление в системе.
Рассчитывается при помощи теоремы вириалов:
$pV=Nk_bT + \frac{\large 1}{\large 3} \left\langle\sum_{i=1}^N \mathbf{r}_i \cdot \mathbf{F}_i \right\rangle \,,$
где $V$ — объём моделируемой области, а угловые скобки обозначают осреднение по времени.

Начальные условия.
Для моделирования канонического ансамбля частиц начальное положение моделируемых молекул может быть задано достаточно произвольно. Для упрощения начальные положения молекул модно задать как центры одинаковы кубических ячеек, заполняющих расчётную область.

Для удовлетворения заданной температуры, компоненты начальных скоростей молекул проще всего задать случайным образом из диапазона [-1;1], затем произвести сдвиг, чтобы соблюсти условие покоя центра масс:
$\sum_{i=1}N \mathbf{v}_i = 0$
После этого скорости масштабируются с учётом заданной температуры.

Виды граничных условий.

Периодическое.
Данное граничное условие позволяет организовать моделирование псевдобесконечной области.
Условие реализуется, если молекула, покидающую расчётную область, появляется с противоположной стороны с сохранением скорости и направления движения.
Термостат.
Данное условие позволяет заведомо поддерживать постоянные температурно-энергетические условия в области.
Условие реализуется таким образом, что достигающая стенки молекула отражается под правильным углом, но скорость её меняется на скорость. определяемую заданной температурой.

Работа в безразмерных переменных
Физические величины, участвующие в молекулярно-динамических расчётах имеют значительные по модулю показатели степени, а потому чрезвычайно неудобны в вычислительном смысле. По этой причине удобно производить обезразмеривание и сводить вычисления к работе безразмерными величинами, при необходимости восстанавливая размерные значения.

Для моделирования поведения аргона удобными оказываются параметры $\sigma=3.405\ \AA$ , $\varepsilon=1.65\cdot10^{-21}$ Дж и $m_0=66.4\cdot10^{-21}$ кг как единицы измерения длины, энергии и массы.

Ниже перечисляется переход к безразмерным параметрам, записываемым со звездочкой, от размерных с использованием выбранной триады основных единиц:
длина $L^*=L/\sigma$ ;
масса $m^*=m/m_0$ ;
энергия $E^*=E/\varepsilon$ ;
объем $V^*=V/\sigma^3$ ;
плотность $\rho^*=\rho\sigma^3/m_0$ ;
скорость $v^*=v\sqrt{m_0/\varepsilon}$ ;
ускорение $a^*=a\sigma m_0/\varepsilon$ ;
сила $F^*=F\sigma/\varepsilon$ ;
температура $T^*=Tk_b/\varepsilon$ ;
давление $p^*=p\sigma^3/\varepsilon$ ;
время $t^*=t\sqrt{\varepsilon/m_0\sigma^2}$ .

Молекулярная динамика

Навигация