Ich versuche herauszufinden, was ein solches Problem verursacht haben könnte: Ich simuliere eine binäre Lennard-Jones-Mischung (zwei Arten von Atomen, die über das LJ-Potential interagieren) im NVE-Ensemble mit dem Geschwindigkeits-Verlet-Algorithmus. Das Problem, das ich habe, ist, dass meine kinetische Energie allmählich (keine scharfe Änderung bei jedem einzelnen Schritt) auf Null geht und die potenzielle Energie allmählich (keine scharfe Änderung bei jedem einzelnen Schritt) auf einen negativen Wert abfällt, was offensichtlich nicht korrekt ist . Ich habe die Positionen und Geschwindigkeiten der Atome ausgedruckt, und ich finde, dass die Geschwindigkeiten zu verschwindenden Werten gehen und sich die Positionen fast nicht ändern, wenn sie sehr lange laufen. Die Randbedingung ist gut, also ist die Lautstärke noch fest. Das bedeutet im Grunde, dass die Atome in meiner Simulation irgendwann zum Stillstand kommen.
Es ist leicht zu überprüfen, ob mein Geschwindigkeits-Verlet-Algorithmus korrekt ist, aber was könnte Ihrer Meinung nach die Ursache für diese Art von seltsamen Ergebnissen sein? Es sieht so aus, als würde sich die Energie in meiner Simulation auflösen. Der Code zur Berechnung der potentiellen Energie und der Kraft ist hier, mit einem Cutoff beim Abstand rc, wobei die Parameter die Kob-Andersen-Parameter sind. Nklein ist die Anzahl der kleineren Atome, während N die Gesamtzahl der Atome ist. Für jedes Array (r, f) sind die ersten (N-Nkleinen) Elemente die Werte für die größeren Atome, die restlichen Nkleinen Elemente für die kleineren Atome.
double total_e ( double * rx, double * ry, double * rz,
double * fx, double * fy, double * fz,
int N, int Nsmall, double L) {
int i,j;
double dx, dy, dz, r2, rc2, r6i, rc6i, epsilon, sigma;
double e = 0.0, hL=L/2.0,f;
/* Zero the forces */
for (i=0;i<N;i++) {
fx[i]=fy[i]=fz[i]=0.0;
}
for (i=0;i<(N-1);i++) {
for (j=i+1;j<N;j++) {
/* Kob Andersen parameters here */
if (i < N - Nsmall && j < N - Nsmall){
epsilon = 1.0;
sigma = 1.0;
} else if (i >= N - Nsmall && j >= N - Nsmall){
epsilon = 0.5;
sigma = 0.88;
} else {
epsilon = 1.5;
sigma = 0.8;
}
dx = (rx[i]-rx[j]);
dy = (ry[i]-ry[j]);
dz = (rz[i]-rz[j]);
/* Periodic boundary conditions: Apply the minimum image
convention; note that this is *not* used to truncate the
potential as long as there an explicit cutoff. */
if (dx>hL) dx-=L;
else if (dx<-hL) dx+=L;
if (dy>hL) dy-=L;
else if (dy<-hL) dy+=L;
if (dz>hL) dz-=L;
else if (dz<-hL) dz+=L;
r2 = dx*dx + dy*dy + dz*dz;
rc2 = 6.25*sigma*sigma; /* Kob Andersen parameter for rcut*/
if (r2<rc2) {
r6i = 1.0/(r2*r2*r2);
rc6i = 1.0/(rc2*rc2*rc2); // cutoff correction to the potential and to the force
e += 4*epsilon*(pow(sigma,12)*r6i*r6i - pow(sigma,6)*r6i) - 4*epsilon*(pow(sigma,12)*rc6i*rc6i - pow(sigma,6)*rc6i) + 48*epsilon*(pow(sigma,12)*rc6i*rc6i/(sqrt(rc2))-0.5*pow(sigma,6)*rc6i/(sqrt(rc2)))*(sqrt(r2) - sqrt(rc2));
f = 48*epsilon*(pow(sigma,12)*r6i*r6i-0.5*pow(sigma,6)*r6i) - 48*epsilon*(pow(sigma,12)*rc6i*rc6i-0.5*pow(sigma,6)*rc6i);
/* Newton's third law*/
fx[i] += dx*f/r2;
fx[j] -= dx*f/r2;
fy[i] += dy*f/r2;
fy[j] -= dy*f/r2;
fz[i] += dz*f/r2;
fz[j] -= dz*f/r2;
}
}
}
return e;
}
Übrigens, eine andere Frage ist, sollte die Gesamtenergie in einer numerischen Geschwindigkeits-Verlet-Simulation streng konstant sein?
fx, fy, fzEnergie wird in jeder Dynamiksimulation nicht unbedingt genau eingespart, da der numerische Algorithmus den Wert einer Größe (in Ihrem Code anscheinend die Kräfte ) während jedes Zeitschritts "abtastet" und davon ausgeht, dass er über die gesamte endliche Länge konstant bleibt Schritt. In der realen Situation ändert sich die Kraft während des Zeitschritts kontinuierlich.
Ich habe keine Erfahrung mit Molekulardynamik-Simulationen, aber diese Wikiepdia-Seite enthält ein paar Verweise auf das Problem der Energiedrift: https://en.wikipedia.org/wiki/Energy_drift
Die Wiki-Seite gibt eine andere (aber im Wesentlichen äquivalente) Erklärung dieses Verhaltens aus dem ersten Absatz dieser Antwort: Der Verlet-Algorithmus spart zwar Energie, aber die diskretisierten Gleichungen, die Sie lösen, entsprechen einem anderen Hamiltonian als in der realen Situation. Im Allgemeinen hängt der Unterschied zwischen den beiden Hamiltonoperatoren von der Zeit und der Entwicklung der numerischen Lösung ab. Obwohl es für jeden einzeln betrachteten Zeitschritt "gegen Null konvergiert", da die Zeitschrittgröße gegen Null tendiert, konvergiert das Integral des Fehlers über die gesamte Lösung möglicherweise nicht gegen Null oder konvergiert langsamer als eine Fehleranalyse von einem einzigen Schritt vorschlagen.
„Energieeinsparung“ ist womöglich nicht alles. Zum Beispiel gibt es wohlbekannte numerische Integrationsverfahren, die Energie für einfache Probleme wie einen linearen harmonischen Oszillator "exakt" sparen, aber der konstante Wert der Energie ist eine Funktion der Zeitschrittgröße und der falschen Energiemenge in der numerischen System ergibt einen Fehler in der Schwingungsfrequenz gegenüber der ursprünglichen Differentialgleichung - dh die Frequenz der numerischen Lösung ist nicht exakt sondern eher Wo ist der numerische Zeitschritt.
QMechaniker
dmckee --- Ex-Moderator-Kätzchen
Damien