Метод энтропического моделирования

Задачи равновесной статистической термодинамики классических систем можно свести к вычислению статистического интеграла. Например, в каноническом ансамбле:

${\displaystyle Z(N,V,T)={\frac {1}{N!(2\pi \hbar )^{3N}}}\int \exp\{-\beta E(p,q)\}dpdq,\qquad {\mbox{где }}\beta ={\frac {1}{kT}},}$

${\displaystyle N}$ - число частиц, находящихся в объёме ${\displaystyle V}$ при температуре ${\displaystyle T}$, ${\displaystyle \beta =1/kT}$; ${\displaystyle E(p,q)}$ - полная механическая энергия частиц; ${\displaystyle p,q}$ - набор их импульсов и координат, при чем ${\displaystyle dp=d^{3}p_{1}\ldots d^{3}p_{N}\,dq=d^{3}q_{1}\ldots d^{3}q_{N}}$. Классическая энергия ${\displaystyle E(p,q)}$ всегда может быть представлена в виде суммы кинетической ${\displaystyle K(p)}$ и потенциальной ${\displaystyle U(q)}$ энергий. Кинетическая энергия есть квадратичная функция от импульсов, и интегрирование по ним может быть произведено в общем виде. В результате получаем:

${\displaystyle Z(N,V,T)={\frac {1}{N!\Lambda ^{3N}}}\int \exp\{-\beta U(q)\}dq}$

где ${\displaystyle \Lambda ={\sqrt {\frac {2\pi \hbar ^{2}}{mkT}}}}$ - тепловая длина волны де Бройля частиц массы ${\displaystyle m}$ при температуре ${\displaystyle T}$. Таким образом, задача сводится к вычислению конфигурационного интеграла

${\displaystyle Q(N,V,T)=\int \exp\{-\beta U(q)\}dq}$

От интегрирования по координатам можно перейти к интегрированию по энергии:

${\displaystyle Q(\beta )=\int \Omega (E)\exp(-\beta E)dE}$

${\displaystyle \Omega (E)=\int \delta (U(q)-E)dq}$

где ${\displaystyle \Omega (E)}$ - объём части конфигурационного пространства, в которой энергия системы лежит в пределах от ${\displaystyle E}$ до ${\displaystyle E+dE}$, ${\displaystyle \delta }$ - дельта функция.

Вычисления по приведённым формулам мы будем выполнять численными методами. Поэтому от интегралов перейдем к интегральным суммам. Диапазон изменения энергии системы ${\displaystyle E_{min}\leq E\leq E_{max}}$ разбивается на конечное число ${\displaystyle (N_{b})}$ равных отрезков. Определяются значения ${\displaystyle \Omega (E_{i}),i=1,2,...,N_{b}}$. В итоге, для любой величины ${\displaystyle f}$ её средние канонические могут быть вычислены по формуле:

${\displaystyle <f>(\beta )={\frac {\sum \limits _{i=1}^{N_{b}}f_{i}\Omega _{i}\exp(-\beta E_{i})}{\sum \limits _{i=1}^{N_{b}}\Omega _{i}\exp(-\beta E_{i})}}}$,

где ${\displaystyle f_{i}}$ — значение величины ${\displaystyle f}$ для ${\displaystyle i}$-го отрезка энергии. Поскольку ${\displaystyle \Omega (E)}$ входит линейно и в числитель, и в знаменатель формулы для ${\displaystyle <f>}$, то ${\displaystyle \Omega (E)}$ можно понимать не только как объём, но и как долю конфигурационного пространства, соответствующую энергии ${\displaystyle E}$. В каждом состоянии (конфигурации) система обладает определённой энергией. Т.е. каждому состоянию (конфигурации) системы можно сопоставить точку на энергетической шкале (оси) в пространстве энергий (это пространство одномерно). Последовательности случайных изменений конфигурации системы соответствует случайное блуждание точки в пространстве энергий. Моделируя процесс случайных блужданий с помощью метода Монте-Карло и зная или вычисляя значения ${\displaystyle \Omega _{i}}$ , мы можем находить средние значения физических величин.

Алгоритм энтропического моделирования основан на следующем обстоятельстве. Совершая случайное блуждание в пространстве энергий с вероятностями перехода, пропорциональными обратной плотности состояний ${\displaystyle 1/\Omega (E)}$, мы получаем равномерное распределение по энергиям. Иными словами, подобрав вероятности перехода такими, что посещение всех энергетических состояний стало бы равномерным, можно получить изначально неизвестную плотность состояний ${\displaystyle \Omega (E)}$.

Напишем конфигурационный интеграл в каноническом ансамбле в виде:

${\displaystyle Q(\beta )=\int e^{-\beta E(q)}dq=\int \Omega (E)e^{-\beta E}dE=\int e^{S(E)-\beta E}dE}$

где ${\displaystyle S(E)=k\ln \Omega (E)}$ - энтропия при заданном значении ${\displaystyle E}$ (иногда ${\displaystyle k}$ будет опускаться, т.к. в моделировании не обязательно учитывать эту константу).

Осуществляя блуждание в конфигурационном пространстве с вероятностями перехода, удовлетворяющими соотношению детального баланса

${\displaystyle {\frac {p(q_{1}\rightarrow q_{2})}{p(q_{2}\rightarrow q_{1})}}=e^{-\beta (E(q_{2})-E(q_{1}))}}$,

получают каноническую выборку состояний ${\displaystyle P(q)\sim e^{-\beta E(q)}}$ (или ${\displaystyle P(E)\sim e^{S(E)-\beta E}}$). Произвольной выборке энергетических состояний ${\displaystyle P(E)\sim e^{A(E)}=e^{S(E)-J(E)}}$, где ${\displaystyle A(E)}$ — произвольная функция, ${\displaystyle J(E)=S(E)-A(E)}$, соответствует условие

${\displaystyle {\frac {p(q_{1}\rightarrow q_{2})}{p(q_{2}\rightarrow q_{1})}}=e^{-\beta (J(E(q_{2}))-J(E(q_{1})))}}$.

При ${\displaystyle J(E)=S(E)}$, в процессе блуждания должна получиться равномерная, в пределах статистического разброса, выборка энергетических состояний, ${\displaystyle P(E)\sim const}$. В этом случае из определения энтропии следует

${\displaystyle {\frac {p(q_{1}\rightarrow q_{2})}{p(q_{2}\rightarrow q_{1})}}={\frac {\Omega (E(q_{1}))}{\Omega (E(q_{2}))}}}$

Таким образом, если при некотором выборе вероятностей перехода получить равномерное посещение энергетических состояний, то можно вычислить плотность состояний ${\displaystyle \Omega (E)}$, а следовательно, и конфигурационный интеграл ${\displaystyle Q(\beta )}$.