Estoy tratando de averiguar qué podría haber causado tal problema: estoy simulando una mezcla binaria de Lennard-Jones (dos tipos de átomos que interactúan a través del potencial LJ) en el conjunto NVE con el algoritmo Verlet de velocidad. El problema que tengo es que mi energía cinética va gradualmente (sin cambios bruscos en ningún paso individual) a cero, y la energía potencial cae gradualmente (sin cambios bruscos en ningún paso individual) a un valor negativo, lo que obviamente no es correcto. . Imprimí las posiciones y velocidades de los átomos, y descubrí que las velocidades van a valores que se desvanecen y las posiciones casi no cambian si dura mucho tiempo. La condición de contorno es buena, por lo que el volumen sigue siendo fijo. Básicamente significa que en mi simulación, los átomos finalmente se detienen.
Es fácil verificar que mi algoritmo Verlet de velocidad es correcto, pero ¿cuál crees que podría ser la causa de este tipo de resultados extraños? Parece que la energía se está disipando en mi simulación. El código para calcular la energía potencial y la fuerza está aquí, con un corte en la distancia rc, siendo los parámetros los parámetros de Kob Andersen. Nsmall es el número de átomos más pequeños, mientras que N es el número total de átomos. Para cada arreglo (r, f), los primeros elementos (N-Npequeños) son los valores para los átomos más grandes, con el resto Npequeños elementos para los átomos más pequeños.
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;
}
Por cierto, otra pregunta es, ¿debería la energía total ser estrictamente constante en una simulación numérica Verlet de velocidad?
La energía no necesariamente se conservará exactamente en cualquier simulación dinámica, porque el algoritmo numérico está "muestreando" el valor de alguna cantidad (en su código, aparentemente las fuerzas fx, fy, fz) durante cada paso de tiempo y suponiendo que permanece constante en toda la longitud finita paso. En la situación del mundo real, la fuerza cambia continuamente durante el paso de tiempo.
No tengo experiencia con simulaciones de dinámica molecular, pero esta página de Wikiepdia tiene un par de referencias al problema de la deriva de energía: https://en.wikipedia.org/wiki/Energy_drift
La página Wiki brinda una explicación diferente (pero esencialmente equivalente) de este comportamiento del primer párrafo de esta respuesta: el algoritmo de Verlet conserva energía, pero las ecuaciones discretizadas que está resolviendo corresponden a un hamiltoniano diferente de la situación del mundo real. En general, la diferencia entre los dos hamiltonianos dependerá del tiempo y de la evolución de la solución numérica. A pesar de que "converge a cero" para cada paso de tiempo considerado individualmente, como el tamaño del paso de tiempo tiende a cero, la integral del error en toda la solución puede no converger a cero, o convergerá a un ritmo más lento que un análisis de error de un solo paso sugieren.
La "conservación de la energía" puede no ser la historia completa en cualquier caso. Por ejemplo, existen métodos de integración numérica bien conocidos que conservan la energía "exactamente" para problemas simples como un oscilador armónico lineal, pero el valor constante de la energía es una función del tamaño del paso de tiempo y la cantidad incorrecta de energía en el valor numérico. El sistema da como resultado un error en la frecuencia de oscilación en comparación con la ecuación diferencial original, es decir, la frecuencia de la solución numérica no es exactamente sino más bien dónde es el paso de tiempo numérico.
qmecanico
dmckee --- gatito ex-moderador
Damián