Как определить равновесие в моделировании Монте-Карло NVTNVTNVT?

Я вижу, что этот вопрос снова был поднят на главную страницу. Недавно я ответил на аналогичный вопрос на https://math.stackexchange.com/a/2920136/575517 , поэтому я изложу его суть здесь, если кто-то сочтет это полезным.

Этот вопрос «Как я узнаю, что прогон симуляции достиг равновесия?» часто приукрашивают и оставляют как «эмпирическое правило». Как обсуждалось в других ответах, обычно требуется, чтобы время «уравновешивания» или «прижигания» было не менее продолжительным, чем время корреляции.$\tau_A$переменной$A$представляющие интерес или, что еще лучше, все интересующие переменные (занимающие самые длинные$\tau$). Эту часть траектории отбрасывают и после этого начинают накапливать средние значения.

Проблема здесь в том, что нужна оценка$\tau_A$заранее, а также что этот аргумент в общих чертах основан на теории линейного отклика : а именно, что релаксация слегка возмущенного состояния к равновесию происходит в масштабе времени, определяемом формулой$\tau_A$, которое является свойством корреляционной функции равновесного времени. Нет никакой гарантии, что релаксация из произвольно подготовленной начальной конфигурации будет следовать этому закону, даже если$\tau_A$может дать разумное руководство.

Однако мне известно по крайней мере об одной статье, в которой была предпринята попытка Джона Чодеры объективно подойти к этому вопросу: https://doi.org/10.1101/021659 , которая также была опубликована в J Chem Theo Comp, 12, 1799 . (2016) .

Я не буду пытаться воспроизвести здесь математику, но основная идея состоит в том, чтобы использовать процедуру оценки статистических ошибок в коррелированных последовательностях данных, которая включает оценку времени корреляции (или статистической неэффективности, которая представляет собой интервал между эффективно независимыми выборками). ) - и применяя его к интервалу$(t_0,t_{\text{max}})$который охватывает период между (предлагаемым) окончанием периода уравновешивания,$t_0$, и конец всего набора данных,$t_{\text{max}}$. Этот расчет ведет себя предсказуемым образом, если набор данных находится в равновесии: колебания$\langle\delta A_{\Delta t}^2\rangle$конечного среднего

$$\delta A_{\Delta t} = A_{\Delta t}-\langle A\rangle, \qquad A_{\Delta t} = \frac{1}{\Delta t} \int_t^{t+\Delta t} dt \, A(t)$$ известным образом зависят от интервала усреднения$\Delta t$. Истоки времени$t$выбираются в пределах интересующего интервала,$(t_0,t_{\text{max}}-\Delta t)$; колебания рассчитываются по всем периодам$\Delta t$лежащие в этом промежутке. Метод систематически выполняет этот расчет в зависимости от$t_0$, уменьшая его от начального высокого значения к началу прогона. В какой-то момент ожидается, что будут замечены отклонения от ожидаемого поведения. Этот момент считается концом периода равновесия.

В любом случае, чтение этого документа должно помочь в прояснении этого вопроса. Автор также предоставляет часть программного обеспечения Python для автоматического выполнения вычислений, поэтому оно также может быть полезным на практике.

Как говорится во многих комментариях, единого и лучшего ответа не существует, каждый использует свой метод. Решение, которое вы нашли, хорошее, но как определить, что равновесие достигнуто?

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

$$N = 10$$

И с параметрами, которые определяют ваше равновесие, вы вычисляете среднее значение и стандартное отклонение:

$$\langle P\rangle, \langle \Delta P^2\rangle$$

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

$$Var(\langle P\rangle_i , \langle P\rangle_{i-1},...,\langle P\rangle_{i-n})\rightarrow0$$

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

П.Д.: Флуктуации после того, как система достигла равновесия, нормальны, а также эти флуктуации увеличиваются с повышением температуры.

Вы не можете знать наверняка.

При моделировании система постоянно хаотично изменяется, и в целом возможно, что после некоторого периода относительного постоянства макроскопических переменных произойдет какое-то радикальное изменение, такое как фазовый переход.

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

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

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

Я считаю совершенно нормальным иметь большие колебания термодинамических функций при повышении температуры. В чем я не уверен, так это в том, относится ли это к реальной ситуации или нет...

Если вы знаете форму результирующей функции распределения частиц, вы можете догадаться, достигла ли ваша система равновесия путем вычисления ее среднего значения, производя расчет функции распределения через каждые n временных шагов (300 временных шагов было хорошо в моих симуляциях Монте-Карло).

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

Если вы настраиваете систему для большого диапазона температур, изменение параметра смещения может быть хорошей идеей, поскольку скорость принятия зависит от отношения (новая энергия к старой энергии)/T.

Я нашел ответ в книге Моделирование методом Монте-Карло в статистической физике - введение (К.Биндер, Д.В.Херманн), стр. 35.

Чтобы определить равновесие, нам нужно запустить симуляцию несколько раз, скажем,$n_{run}$раз. мы определяем среднее$<>_T$в среднем после$t$этапы моделирования:

$$\langle A(t)\rangle_T = \frac{1}{n_{run}} \sum_{l=1}^{n_{run}}{A(t,l)}$$ Где А — какое-то физическое свойство, например энергия или давление.$A(t,n)$- значение этого свойства во времени (шаге моделирования) t n-го моделирования. Теперь определим «нелинейную функцию релаксации»: $$\phi_A(k)=\frac{\langle A(k)\rangle_T - \langle A(\infty)\rangle_T}{\langle A(0)\rangle_T - \langle A(\infty)\rangle_T}$$ Где$\langle A(\infty)\rangle_T$является средним значением A на последнем шаге моделирования, для$n_{run}$симуляции. Эта функция релаксации обращается в нуль после некоторого шага t. время релаксации этой функции равно: $$\tau_A=\int_0^{\infty}\phi(t)dt$$ до сих пор не уверен, почему это время релаксации, я был бы рад, если бы кто-нибудь мог мне это объяснить. Итак, если это время релаксации, то время, за которое система достигает равновесия, должно быть намного больше, чем это время: $$t_{simulation} \gg \tau_A$$ поскольку разные свойства могут иметь разное время релаксации, мы должны учитывать самое медленное время релаксации, чтобы привести систему в равновесие.