Независимость от временного шага кода молекулярной динамики

Я пишу базовый код MD на C ++, используя потенциал LJ для системы NVE. Начальная конфигурация — FCC, а начальные скорости генерируются случайным образом.

Я сталкиваюсь со странной проблемой, заключающейся в том, что эволюция системы, по-видимому, не зависит от временного шага, который я выполняю. Насколько я понимаю, потери энергии меньше для маленьких временных шагов и больше для больших временных шагов. Однако я получаю тот же результат в конце моделирования с точки зрения энергии, бегу ли я (0,0001 шага) * (10000 шагов) или 0,001 * 1000 и так далее.

Весь код слишком большой для того, чтобы я мог публиковать его здесь, поэтому я публикую то, что считаю уместным, и опускаю биннинг и т. Д., Пожалуйста, дайте мне знать, если потребуется какая-либо дополнительная информация. Я прошел через бесчисленное множество кодов, доступных онлайн, и хотя они выглядят похожими на мои, я просто не могу понять, в чем разница и где я ошибаюсь.

Основной cpp содержит следующий цикл

for (int i=0; i<t;i++)
{
potential_calc(neighlist,fromfile, run_parameters,i);//calculating the force fields
velverlet(neighlist,fromfile, run_parameters, bin, dt);//calculating the velocities
}

Объявления 2-х cpp файлов для потенциального расчета & интеграция Verlet

void potential_calc(neighborlist_type *neighlist, config_type *fromfile, potential *run_parameters, int t)

void velverlet(neighborlist_type *neighlist, config_type *fromfile, potential *run_parameters, bin_type *bin, double dt)

Код для расчета силы — потенциал_кала.cpp ниже

for (long i=0; i<fromfile->N; i++)
{
long atom_p=i;for (long j=0; j<neighlist[i].countsn; j++)
{

long atom_s=neighlist[i].numb[j];

for (int k=0; k<Dim; k++)
{
dist[k]= fromfile->r[atom_p][k] - (fromfile->r[atom_s][k] + neighlist[atom_p].xyz[j][k]*fromfile->L[k]);
//the .xyz indicates the image being considered real or mirror(if mirror then in which direction)
}
disp2 = pow(dist[0],2)+pow(dist[1],2)+pow(dist[2],2);
if (disp2<rb2)
{
int c1=fromfile->c[atom_p];
int c2=fromfile->c[atom_s];
double long force_temp;
disp=pow(disp2,0.5);

sig_r6=pow(run_parameters->sigma[c1-1][c2-1]/disp,6);//(sigma/r)^6
sig_r8=pow(run_parameters->sigma[c1-1][c2-1]/disp,8);//(sigma/r)^8

run_parameters->pe[atom_p] += (4*run_parameters->epsilon[c1-1][c2-1]*((sig_r6*sig_r6)-sig_r6)) - potential_correction[c1-1][c2-1];

force_temp=(-1*((48*run_parameters->epsilon[c1-1][c2-1])/pow(run_parameters->sigma[c1-1][c2-1],2)*((sig_r6*sig_r8)-((sig_r8)*0.5))));
for (int k=0; k<Dim;k++)
{
run_parameters->force[atom_p][k]+=force_temp*(-1*dist[k]);
}
}
}
//calculating kinetic energy
run_parameters->ke[atom_p] = 0.5*(pow(fromfile->v[atom_p][0],2)+pow(fromfile->v[atom_p][1],2)+pow(fromfile->v[atom_p][2],2));

}

Как только вычисление силы выполнено, оно переходит к обновлению скорости и положения в velverlet.cpp

for (long i=0; i<fromfile->N; i++)
{
for (int j=0; j<Dim; j++)
{
fromfile->v[i][j] += (dt*run_parameters->force[i][j]);
}
}

for (long i=0; i<fromfile->N; i++)
{
for (int j=0; j<Dim; j++)
{
fromfile->r[i][j] += dt*fromfile->v[i][j];
}
}

Могут быть небольшие различия в том, как верлет скорости реализован разными людьми, но я не могу понять, как я получаю независимые от шага по времени результаты.

Пожалуйста помоги. Любой вклад приветствуется
Извините, если какое-либо форматирование / пометка неправильные, я впервые публикую здесь