Кинетическая энергия постепенно стремится к нулю в моей скоростной симуляции Verlet MD.

Я пытаюсь выяснить, что могло вызвать такую ​​​​проблему: я моделирую бинарную смесь Леннарда-Джонса (два типа атомов, взаимодействующих через LJ-потенциал) в ансамбле NVE с помощью алгоритма скорости Верле. Проблема, которая у меня есть, заключается в том, что моя кинетическая энергия постепенно (без резких изменений на любом отдельном шаге) стремится к нулю, а потенциальная энергия постепенно падает (без резких изменений на любом отдельном шаге) до отрицательного значения, что, очевидно, неверно. . Я распечатал положения и скорости атомов и обнаружил, что скорости стремятся к нулю, а положения почти не меняются, если работать очень долго. Граничное условие хорошее, поэтому объем по-прежнему фиксирован. В основном это означает, что в моей симуляции атомы в конце концов останавливаются.

Легко проверить правильность моего алгоритма скорости Верле, но как вы думаете, что может быть причиной таких странных результатов? Похоже, что в моей симуляции энергия рассеивается. Код для расчета потенциальной энергии и силы находится здесь, с отсечкой на расстоянии rc, параметры являются параметрами Коба Андерсена. Nsmall — это количество более мелких атомов, а N — это общее количество атомов. Для каждого массива (r, f) первые (N-Nsmall) элементы являются значениями для более крупных атомов, а остальные Nsmall-элементы — для меньших атомов.

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;

}

Кстати, другой вопрос: должна ли полная энергия быть строго постоянной при численном моделировании скорости Верле?