Доклады Российской академии наук. Физика, технические науки, 2020, T. 494, № 1, стр. 5-9

АППРОКСИМАЦИЯ ТРАЕКТОРИИ В ФАЗОВОМ ПРОСТРАНСТВЕ МЕТОДОМ ЭФФЕКТИВНЫХ МОД И ОЦЕНКА ТЕМПЕРАТУРЫ

По аналогии с классической термодинамикой, температура конечномерных систем, какими являются и кластеры, оценивается по следующей формуле [3]:

где E – средняя по траектории кинетическая энергии кластера, k_B – константа Больцмана, N_df  – число внутренних степеней свободы. Для описания динамики кластера фазовую траекторию системы удобно представить как суперпозицию независимых (в некотором смысле) составляющих, каждая из которых описывает некоторое коллективное движение его атомов. Например, для системы из n частиц в приближении малых колебаний вблизи состояния равновесия фазовую траекторию можно представить в виде линейной комбинации нормальных мод. Динамика мод описывается системой независимых колебательных уравнений [6]. Каждая нормальная мода отражает коллективное колебание системы материальных точек с фиксированной частотой. В случае невращающегося кластера в системе центра масс число нормальных мод N_nm оценивается следующим образом: N_nm = 3n – Q – 6 (здесь n – число атомов, Q – число жестких связей в кластере). Для октамера воды с PIT4P потенциалом N_nm = 42 (в каждой молекуле воды зафиксированы 3 внутренние координаты: расстояния между атомами водорода и кислорода, а также угол между линиями ковалентных связей в молекуле).

При росте полной энергии системы линеаризованная модель становится слишком грубым приближением. В этом случае мы используем подход, описанный в работе [7], представляющий собой модификацию метода главных компонент и состоящий в следующем.

Пусть $x(t) = (q(t),p(t)) \in {{R}^{{6n}}}$, $t \in [0,\tau ]$, – фазовая траектория кластера из n частиц. Метод главных компонент позволяет представить траекторию в виде разложения x(t) = $\sum\limits_{i = 1}^{6n} {(x(t),{{e}_{i}})} \,{{e}_{i}}$ по ортонормированному базису $\left\{ {{{e}_{i}}} \right\} \subset {{R}^{{6n}}}$, обладающему следующим экстремальным свойством: любая m-я частичная сумма

приближает траекторию x(t) в среднем по времени не хуже, чем m-я частичная сумма $x_{m}^{'}(t)$ = $\sum\limits_{i = 1}^m {(x(t),e_{i}^{'})\,e_{i}^{'}} $ в разложении по любому другому ортонормированному базису $\{ e_{i}^{'}\} \subset {{R}^{{6n}}}$. Векторы базиса {e_i} есть собственные векторы линейного оператора X : R⁶ⁿ → → ${{R}^{{6n}}}$, определенного для любого $y \in {{R}^{{6n}}}$ равенством

Они соответствуют собственным числам оператора X, упорядоченным по невозрастанию: $X{{e}_{i}} = \alpha _{i}^{2}{{e}_{i}}$, i = 1, …, 6n, $\alpha _{1}^{2} \geqslant \alpha _{2}^{2} \geqslant \; \ldots \; \geqslant \alpha _{{6n}}^{2} \geqslant 0$. При этом

т.е. m-я частичная сумма (2) в среднем по времени приближает траекторию x(t) по норме с минимальной погрешностью $\sum\limits_{i = m + 1}^{6n} {\alpha _{i}^{2}} $, и доля среднего по времени квадрата нормы траектории x(t), содержащейся в m-й частичной сумме (2), максимальна и равна

Таким образом, траектория x(t), $t \in [0,\tau ]$, может быть представлена в виде суммы эффективных мод $(x(t),{{e}_{i}})\,{{e}_{i}}$, векторы ${{e}_{1}},{{e}_{2}},\; \ldots ,\;{{e}_{{6n}}}$ представляют собой взаимно ортогональные в R⁶ⁿ коллективные моды движения кластера, а зависящие от времени скалярные произведения $(x(t),{{e}_{i}})$ – амплитуды эффективных мод. Средний квадрат i-й эффективной моды равен $\alpha _{i}^{2}$, i = 1, …, 6n, а его доля в величине среднего по времени квадрата нормы траектории x(t) равна $\frac{{\alpha _{i}^{2}}}{{\sum\limits_{i = 1}^{6n} {\alpha _{i}^{2}} }}$.

В настоящем сообщении вместо траектории $x(t) \in {{R}^{{6n}}}$, $t \in [0,\tau ]$, предлагается рассматривать функцию от фазовых координат ${{E}^{{1/2}}}(t) \in {{R}^{{3n}}}$, квадрат нормы которой в каждый момент времени равен кинетической энергии кластера. Если движение молекулярного кластера происходит в трехмерном декартовом пространстве, то эта функция получена путем деления каждой координаты импульса каждой частицы на квадратный корень из ее удвоенной массы, так что ее координаты запишутся в виде $\left( {\frac{{{{p}_{1}}(t)}}{{\sqrt {2{{m}_{1}}} }}} \right.,\frac{{{{p}_{2}}(t)}}{{\sqrt {2{{m}_{1}}} }},\frac{{{{p}_{3}}(t)}}{{\sqrt {2{{m}_{1}}} }},\frac{{{{p}_{4}}(t)}}{{\sqrt {2{{m}_{2}}} }},\; \ldots ,\;\left. {\frac{{{{p}_{{3n}}}(t)}}{{\sqrt {2{{m}_{n}}} }}} \right)$. Тогда, заменив в формулах (2), (3) $x(t) \in {{R}^{{6n}}}$ на E^1/2(t) ∈ R³ⁿ, решив задачу $E{{g}_{i}} = \lambda _{i}^{2}{{g}_{i}}$ на собственные значения линейного оператора $E:{{R}^{{3n}}} \to {{R}^{{3n}}}$, такого, что

и найдя собственные векторы g_i, соответствующие его собственным числам, упорядоченным согласно неравенствам $\lambda _{1}^{2} \geqslant \lambda _{2}^{2} \geqslant \; \ldots \; \geqslant \lambda _{{3n}}^{2} \geqslant 0$, получим, что траекторию $q(t) \in {{R}^{{3n}}}$ (или $p(t) \in {{R}^{{3n}}}$) можно разложить на ортогональные в R³ⁿ составляющие (моды) q(t) = $\sum\limits_{i = 1}^{3n} {(q(t),{{g}_{i}})\,} {{g}_{i}}$, p(t) = $\sum\limits_{i = 1}^{3n} {(p(t),{{g}_{i}})} \,{{g}_{i}}$, упорядоченные по убыванию максимальной доли средней по времени кинетической энергии, соответствующей этой моде.

Траектория движения кластера как целого рассчитывалась методом молекулярной динамики с начальными условиями, полученными методом Монте-Карло [8]. Начальное состояние системы соответствовало энергии глобального минимума Н₀ = –3.17 эВ и S₄ – конфигурации кластера. Для описания взаимодействия между молекулами воды использовался аналитический потенциал жесткого типа TIP4P [9].

РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

На рис. 1 показано, каким образом средняя по траектории кинетическая энергия E кластера распределяется по эффективным модам коллективного движения частиц системы. Три кривые отражают зависимость этого распределения от начальной полной энергии H кластера. На рис. 1 хорошо виден процесс активации эффективных мод с ростом энергии системы. Время наблюдения за системой составляло 100 пс в каждом из трех случаев.

Имея распределение кинетической энергии кластера по эффективным модам (аналогичное рис. 1), можно получить энтропийную оценку числа эффективных мод по следующей формуле: N = $\exp \left( { - \sum\limits_{j = 1}^{3n} {{{\theta }_{j}} \times \ln ({{\theta }_{j}})} } \right)$. Здесь ${{\theta }_{j}} = \frac{{\left\langle {{{E}_{j}}} \right\rangle }}{E}$ – доля кинетической энергии в j-й моде: $\left\langle {{{E}_{j}}} \right\rangle $ = = $\frac{1}{\tau }\int\limits_0^\tau {{{{({{E}^{{1/2}}}(t),{{g}_{i}})}}^{2}}dt} $, нормированная на $E$. Для квантовых систем асимптотика энтропии представлена в [10]. На рис. 2 показаны результаты расчета числа эффективных мод N кластера с ростом полной энергии H. Видно, что число эффективных мод не превышает число нормальных мод N_nm, которое для октамера воды с жестким потенциалом взаимодействия равно 42, только в области низких энергий кластера.

Резкий скачок числа эффективных мод системы в области H ∼ –2.2 эВ приводит к падению температуры T_eff, рассчитанной по формуле (1) с учетом N (рис. 3). Следует заметить, что аналогичное немонотонное поведение температуры было обнаружено и для других кластеров [11], и является характерным для нелинейных конечномерных систем в области фазового перехода.

На рис. 4 представлены зависимости числа эффективных мод от времени наблюдения за траекторией при разных значениях полной энергии кластера. Так как число эффективных мод октамера воды растет со временем, в отличие от числа нормальных мод, которое остается постоянным, система может рассматриваться как линейная только для малых значений t. C течением времени нелинейность системы становится более значимой, в результате траектория покидает линейное подпространство фазовой траектории, определенное линеаризацией.

ЗАКЛЮЧЕНИЕ

Результаты, представленные в данной работе, показывают, что

1) число эффективных мод коллективной динамики кластера в фазовом пространстве существенно зависит от полной энергии системы, в отличие от числа нормальных мод, которое остается неизменным при любой его энергии;

2) поведение температуры с ростом энергии, полученное двумя методами, позволяет оценить пороговое значение энергии кластера, ниже которого метод нормальных мод может быть использован;

3) зависимости числа эффективных мод кластера от времени наблюдения за траекторией в фазовом пространстве дают временн${\text{ы}}'$е оценки линейной динамики системы;

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