Wie bestimme ich das Gleichgewicht in einer Monte-Carlo-NVTNVTNVT-Simulation?

Wie ich sehe, wurde diese Frage erneut auf die Homepage verschoben. Ich habe kürzlich eine ähnliche Frage unter https://math.stackexchange.com/a/2920136/575517 beantwortet , also werde ich hier die Essenz davon geben, falls es jemand hilfreich findet.

Diese Frage "Woher weiß ich, dass der Simulationslauf das Gleichgewicht erreicht hat?" wird oft beschönigt und als „Faustregel“ belassen. Wie in anderen Antworten besprochen, wird normalerweise verlangt, dass eine "Äquilibrierungs-" oder "Einbrennzeit" mindestens so lang wie die Korrelationszeit sein sollte$\tau_A$der Variablen$A$von Interesse oder, noch besser, alle interessierenden Variablen (dauert am längsten$\tau$). Man verwirft diesen Teil der Trajektorie und beginnt danach, Mittelwerte zu akkumulieren.

Probleme hier sind, dass man eine Schätzung braucht$\tau_A$im Voraus, und auch, dass dieses Argument lose auf der Theorie der linearen Reaktion basiert : nämlich, dass die Entspannung eines leicht gestörten Zustands bis zum Gleichgewicht auf einer durch gegebenen Zeitskala erfolgt$\tau_A$, was eine Eigenschaft der Gleichgewichtszeitkorrelationsfunktion ist. Es gibt jedoch keine Garantie dafür, dass die Entspannung von einer willkürlich erstellten Anfangskonfiguration diesem Gesetz folgt$\tau_A$könnte eine vernünftige Anleitung geben.

Mir ist jedoch mindestens ein Artikel bekannt, in dem ein Versuch von John Chodera unternommen wurde, ihn objektiv anzugehen: https://doi.org/10.1101/021659 , der auch in J Chem Theo Comp, 12, 1799 veröffentlicht wurde (2016) .

Ich werde hier nicht versuchen, die Mathematik zu reproduzieren, aber die Grundidee besteht darin, das Verfahren zum Schätzen statistischer Fehler in korrelierten Datenfolgen zu verwenden - was das Schätzen der Korrelationszeit (oder der statistischen Ineffizienz, dh des Abstands zwischen effektiv unabhängigen Stichproben) beinhaltet ) - und auf das Intervall anwenden$(t_0,t_{\text{max}})$die den Zeitraum zwischen dem (vorgeschlagenen) Ende der Äquilibrierungsperiode abdeckt,$t_0$, und das Ende des gesamten Datensatzes,$t_{\text{max}}$. Diese Berechnung verhält sich auf vorhersagbare Weise, wenn sich der Datensatz im Gleichgewicht befindet: die Schwankungen$\langle\delta A_{\Delta t}^2\rangle$eines endlichen Zeitmittels

$$\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)$$ hängen in bekannter Weise vom Mittelungsintervall ab$\Delta t$. Die Zeit entsteht$t$innerhalb des interessierenden Intervalls gewählt werden,$(t_0,t_{\text{max}}-\Delta t)$; die Schwankungen werden aus allen Perioden berechnet$\Delta t$in diesem Intervall liegen. Das Verfahren führt diese Berechnung systematisch in Abhängigkeit von durch$t_0$, wodurch er von einem anfänglich hohen Wert zu Beginn des Laufs verringert wird. An einem gewissen Punkt wird erwartet, dass Abweichungen vom erwarteten Verhalten zu sehen sind. Dieser Punkt wird als Ende der Äquilibrierungsperiode angesehen.

Wie auch immer, das Lesen dieses Papiers sollte hilfreich sein, um dieses Problem zu klären. Der Autor stellt auch ein Stück Python-Software zur Verfügung, um die Berechnung automatisch zu implementieren, so dass es auch von praktischem Nutzen sein kann.

Wie viele Kommentare sagen, gibt es keine einzige und beste Antwort, jede verwendet eine andere Methode. Die Lösung, die Sie gefunden haben, ist gut, aber wie definieren Sie, wann das Gleichgewicht erreicht ist?

Dazu müssen Sie die letzten Werte der Simulation (Energie, Druck usw.) überprüfen, also wählen Sie eine Reihe früherer Konfigurationen aus, die Sie überprüfen werden:

$$N = 10$$

Und mit den Parametern, die Ihr Gleichgewicht definieren, berechnen Sie den Mittelwert und die Standardabweichung:

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

Diese Werte sollten sich nach einigen Schritten nicht zu sehr ändern. Wenn Sie einige Mittelwerte und ihre Varianz speichern, werden Sie sehen, dass der Mittelwert gegen den Wert des Systems konvergiert und die Varianz über den zeitlichen Mittelwerten bei jedem Schritt gegen Null geht.

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

Sie müssen also als Parameter für das Gleichgewicht wählen, wie viele Schritte Sie für den Mittelwert berücksichtigen und wie viele Mittelwerte Sie verwenden, um die Varianz zu berechnen.

PD: Die Schwankungen, nachdem das System das Gleichgewicht erreicht hat, sind normal und auch, dass die Schwankungen mit der Temperatur zunehmen.

Du kannst es nicht genau wissen.

In der Simulation ändert sich das System ständig chaotisch, und im Allgemeinen ist es möglich, dass nach einer gewissen Zeit relativer Konstanz makroskopischer Variablen einige radikale Änderungen wie ein Phasenübergang auftreten.

Wenn das System jedoch sehr einfach ist, wie z. B. hartes Kugelgas, kann man theoretisch berechnen, wie der Gleichgewichtszustand für gegebene Randbedingungen wie Volumen und Temperatur aussehen sollte. Wenn sich dann die aus der Simulation abgeleiteten Werte diesen theoretischen Werten annähern, können wir sagen, dass wir höchstwahrscheinlich einen Übergang zum Gleichgewicht beobachten und nichts Überraschendes passieren wird.

Wenn wir die Werte von Variablen im thermodynamischen Gleichgewicht nicht kennen, können wir nicht sagen, dass sich das System diesem nähert.

Es gibt jedoch auch die Idee des eingeschränkten thermodynamischen Gleichgewichts, das nur für einen begrenzten Zeitraum definiert ist, an dem wir interessiert sind. Das System befindet sich dann im Gleichgewicht, wenn die Schätzung der durchschnittlichen kinetischen Energie im gesamten Raumbereich des Systems einheitlich ist und beobachtet wird Schwankungen in Zeit und Raum stimmen mit der Boltzmann-Verteilung überein, und dasselbe gilt für andere interessierende physikalische Größen (durchschnittlicher Druck, potenzielle Nettoenergie).

Größere Schwankungen der thermodynamischen Funktionen bei steigender Temperatur halte ich für völlig normal. Wobei ich mir nicht so sicher bin, ob das zu einer realistischen Situation gehört oder nicht...

Wenn Sie die Form der resultierenden Partikelverteilungsfunktion kennen, können Sie erraten, ob Ihr System das Gleichgewicht erreicht hat, indem Sie seinen Mittelwert berechnen, indem Sie die Verteilungsfunktionsberechnung alle n Zeitschritte durchführen (300 Zeitschritte waren a gut n in meinen Monte-Carlo-Simulationen).

Sie müssen berücksichtigen, dass eine Monte-Carlo-Simulation bereits ein Gleichgewichtssystem simulieren soll, sodass Sie bei ihrer Verwendung keine dynamischen Phänomene erwarten sollten (Sie können das Auftreten von Blasen in der Masse einer Mischung nicht beobachten, z Beispiel).

Wenn Sie das System für einen großen Temperaturbereich einstellen, könnte eine Änderung des Verschiebungsparameters eine gute Idee sein, da die Akzeptanzrate vom Verhältnis (neue Energie-alte Energie)/T abhängt.

Eine Antwort fand ich im Buch Monte-Carlo-Simulationen in der statistischen Physik - eine Einführung (von K. Binder, DWHermann), Seite 35.

Um das Gleichgewicht zu bestimmen, müssen wir die Simulation ein paar Mal ausführen, sagen wir mal$n_{run}$mal. Wir definieren den Durchschnitt$<>_T$als Durchschnitt nach$t$Schritte der Simulation:

$$\langle A(t)\rangle_T = \frac{1}{n_{run}} \sum_{l=1}^{n_{run}}{A(t,l)}$$ Wobei A eine physikalische Eigenschaft ist, zum Beispiel die Energie oder der Druck.$A(t,n)$ist der Wert dieser Eigenschaft in Zeit (Simulationsschritt) t der n-ten Simulation. Nun definieren wir die "nichtlineare Relaxationsfunktion": $$\phi_A(k)=\frac{\langle A(k)\rangle_T - \langle A(\infty)\rangle_T}{\langle A(0)\rangle_T - \langle A(\infty)\rangle_T}$$ Wo$\langle A(\infty)\rangle_T$ist der Durchschnitt von A im letzten Simulationsschritt, z$n_{run}$Simulationen. Diese Relaxationsfunktion geht nach einem gewissen Schritt t auf Null. Die Relaxationszeit dieser Funktion ist: $$\tau_A=\int_0^{\infty}\phi(t)dt$$ Ich bin mir immer noch nicht sicher, warum dies die Entspannungszeit ist, ich würde mich freuen, wenn es mir jemand erklären könnte. Wenn dies also die Entspannungszeit ist, dann muss die Zeit, in der das System das Gleichgewicht erreicht, viel größer sein als diese Zeit: $$t_{simulation} \gg \tau_A$$ Da verschiedene Eigenschaften unterschiedliche Relaxationszeiten haben können, müssen wir die langsamste Relaxationszeit berücksichtigen, um das System ins Gleichgewicht zu bringen.