Вариант 2: Реакция полимеризации

Простейшая реакция катализа записывается в виде: \(A+B\xrightarrow{k}A+C\). Динамику концентрации вещества \(A\) и \(B\) при этом можно описать при помощи следующей системы дифференциальных уравнений:

\[\begin{split}\begin{cases} \frac{dB}{dt} = -kAB,\\ \frac{dC}{dt} = kAB. \end{cases}\end{split}\]

В Системе уравнений, \(A=\frac{N_{A}}{V_{0}}\), \(B=\frac{N_{B}}{V_{0}}\), \(C=\frac{N_{C}}{V_{0}}\), \(k\) — кинетическая константа скорости реакции, \(N_{A}\), \(N_{B}\), \(N_{C}\) — количество реагентов типа \(A\), \(B\) и \(C\) соответственно, \(V_{0}\) — объем системы. Решение системы:

\[\begin{split}\begin{cases} B = B_{0}e^{-kAt},\\ C = B_{0}-B=B_{0}(1-e^{-kAt}). \end{cases}\end{split}\]

Реагирующая система может быть описана в виде конечного объема со случайным начальным распределением \(N\) частиц в нем и заданными периодическими граничными условиями. Будем считать, что реакция происходит со скоростью \(p=0.001\) фс \(^{-1}\) когда частицы \(A\) и \(B\) находятся ближе чем \(r_c=4.0\) нм друг к другу. Взаимодействие частиц можно описать при помощи следущего потенциала:

\[V(\vec{r}_i,\vec{r}_j)=\varepsilon\left[С_1\left(\frac{\sigma}{r_{ij}}\right)^{4} + С_2\left(\frac{\sigma}{r_{ij}}\right)^2\right],\]

где \(\varepsilon=1\) кДж/моль — энергетическая постоянная, \(\sigma=1\) нм — равновесное расстояние между частицами, \(С_1\) и \(С_2\) — константы, \(r_{ij}=|\mathbf{r}_i-\mathbf{r}_j|\) — текущее расстояние между частицами. Если задать постоянные \(С_1=1\), \(С_2=-2\), то взаимодействие между частицами будет притягивающим, а при \(С_1=0\), \(С_2=1\) — отталкивающим. Допустим, что продукт реакции \(C\) полимеризуется. Тогда результатом реакции будет полимер, состоящий из молекул типа \(C\). Для описания процесса полимеризации предлагается использовать потенциал c притягивающими постоянными (\(С_1=1\) и \(С_2=-2\)) в случае, если обе частицы имеют тип \(C\) и отталкивающие (\(С_1=0\), \(С_2=1\)) в другом случае.

Движение частиц описывается уравнениями Ланжевена:

\[m_i\frac{d^2\mathbf{r}_i}{dt^2}=-\nabla_iV-\lambda m_i\frac{d\mathbf{r}_i}{dt}+\eta_i(t),\]

где \(m_i = 1.0\) г/моль — масса частиц, \(V({\mathbf{r}_i}) = V_{cov}({\mathbf{r}_i}) + V_{hic}({\mathbf{r}_i})\) — потенциальная энергия системы, \(\nabla_i V\) — градиент потенциала по координате \(i\)-ой частицы. \(\eta_i(t)\) — случайная сила, описывающая соударения с молекулами воды с нормальным распределением:

\[\langle\eta_i(t)\eta_j(t)\rangle = 2\lambda m_i k_BT\delta_{ij}\delta(t-t')\]

Динамику Ланжевена можно рассматривать как динамику Ньютона с добавлением двух дополнительных сил: силы трения, пропорциональной скорости частицы и случайной силы, распределенной по Гауссу. Коэффициент пропорциональности \(\lambda\) называется константой демпфирования и задает силу влияния среды на молекулу. При больших \(\lambda\) система переходит к Брауновскому движению, при малых — к Ньютоновскому.

Для численного решения уравнения Ланжевена можно использовать следующую разностную схему:

\[ \begin{align}\begin{aligned}\begin{split}\begin{split} m_i\frac{\mathbf{v}_{i}^{n+\frac{1}{2}}-\mathbf{v}_{i}^{n-\frac{1}{2}}}{\tau} &= \mathbf{f}_{i}^{n} - \lambda m_i\frac{\mathbf{v}_{i}^{n+\frac{1}{2}}+\mathbf{v}_{i}^{n-\frac{1}{2}}}{2}+\sqrt{\frac{2k_BT\lambda m_i}{\tau}}\mathbf{r}_i^f\\ \frac{\mathbf{r}_{i}^{n+1}-\mathbf{r}_{i}^{n}}{\tau} &= \mathbf{v}_{i}^{n+\frac{1}{2}}\end{split}\\\end{split}\end{aligned}\end{align} \]

Здесь, \(\lambda=50\) пс \(^{-1}\) — коэфициент демпфирования, \(\mathbf{r}_i^f\) — вектор из трех нормально распределенных случайных величин, с дисперсией \(1\) и математическим ожиданием \(0\). Постоянная Больцмана \(k_B=8.31\times10^{-3}\) кДж/K*моль, шаг по времени \(\tau=0.001\) пс, температура \(T=300\) K.

Отображение структуры в программе VMD

VMD не умеет менять тип молекул. Чтобы при визуализации системы цвет атомов менялся при реакции, можно использовать следующий способ. Частицы \(B\) и \(C\) всегда сохраняются в файл координат в максимально возможном их количестве. Например, если в начале симуляции в системе 100 частиц \(B\) и 0 частиц \(C\), в файле координат сохраняется 100 частиц \(B\) и 100 частиц \(C\). При этом для неактивных частиц задаются координаты (0,0,0). В ходе реакции, частицы \(B\) становятся неактивными, а частицы \(C\) — наоборот переходят из неактивных в активные. Координаты же переносятся от прореагирующей частицы к продукту реакции. Таким образом, если отображать только активные частицы, при реакции прореагирующие частицы будут менять тип.

Для визуализации, в строке Selected atoms нужно написать not x = 0, а во вкладке Trajectory выбрать пункт Update Selection Every Frame.