Молекулярная динамика, скорость верлета: расхождение кинетической энергии

Я пытаюсь написать простую MD-программу на C / C ++ (я привык к C, но я пытаюсь изучать C ++, поэтому мой код немного «смешан» … Я знаю, что это неоптимально, и я перейду на полный C ++, как только я это полностью пойму).

Кажется, что все работает, но у меня есть расхождения в кинетической энергии, система не нагревается, и температура (от упора до К) поднимается от порядка (10 ° К) до порядка (10000 ° К) за один шаг.

Я работаю с небольшим временным шагом 0,002 (общее время симуляции: 30), поэтому у меня не должно быть этой огромной ошибки …

Это мой код, если что-то не понятно, я могу попытаться объяснить это лучше

int main(){

…

int n, t, m, i;
double r, K, U, E,P,  totalE, temperature, d, x,y,z, temp;

…

double data[5][PARTICELLE], vel[3][PARTICELLE],dataNew[3][PARTICELLE]; //0,1,2 are x,y,z. 3, 4 for data are Energy and Pressure
double force[3][PARTICELLE], forceNew[3][PARTICELLE];
double velQ[PARTICELLE]; //square velocity

ofstream out(OUTDATA);

//inizio MD

for(t=0; t<PASSI; t++){
//inizialization
K=0;
U=0;
E=0;
P=0;
fill(data[3], data[3]+PARTICELLE, 0); //E=0 for each particle
fill(data[4], data[4]+PARTICELLE, 0);
fill(velQ, velQ+PARTICELLE, 0);
for(i=0; i<3; i++){
fill(force[i], force[i]+PARTICELLE, 0);
fill(forceNew[i], forceNew[i]+PARTICELLE, 0);
}
for(n=0; n<PARTICELLE; n++) { //for on the n_ particle. A step is a move of n=PARTICELLE particles
for (i = 0; i < 3; i++) { //compute vSquare
velQ[n] += vel[i][n] * vel[i][n];
}
K += 0.5 * MASSA * velQ[n]; //compute Kinetic Energy
for(m=0; m<PARTICELLE; m++){ //loop on m!=n to compute F, E, P
if(m!=n){
r=0;
for(i=0; i<3; i++){ //calculation of radius and x,y,z
d = data[i][m] - data[i][n];
d = d - (NINT(d / LATO) * LATO);
if(i==0)x=d;
if(i==1)y=d;
if(i==2)z=d;
r += d * d;
}
//if (t<2) cout << "x y z" << x << " " << y << " " << z << endl;
r=sqrt(r);
if (r < R) {
data[3][n] += energy(r); //update Energy of n
for(i=0; i<3; i++){
if(i==0)temp=x;
if(i==1)temp=y;
if(i==2)temp=z;
force[i][n]+=forza(r,temp); //compute force (cartesian components)
//if(t<2)cout << "force " <<n << " " << m << " "<< force[i][n] << endl;
}
if (m < n)data[4][n] += (-energy(r) * (1 + r)); //pressure
}
}
}
U+=data[3][n]; //total potential energy
P+=data[4][n]; //total pressure
for (i = 0; i < 3; i++) { //Verlet update, part 1
dataNew[i][n] = data[i][n] + vel[i][n] * DeltaT + 0.5 * force[i][n] * DeltaT * DeltaT / MASSA;
}
for(m=0; m<PARTICELLE; m++){ //update force
if(m!=n){
r=0;
for(i=0; i<3; i++){
d = data[i][m] - dataNew[i][n];
d = d - (NINT(d / LATO) * LATO);
if(i==0)x=d;
if(i==1)y=d;
if(i==2)z=d;
r += d * d;
}
r=sqrt(r);
if (r < R) {
for(i=0; i<3; i++) {
if (i == 0)temp = x;
if (i == 1)temp = y;
if (i == 2)temp = z;
forceNew[i][n] += forza(r, temp);

}
}
}
}
for(i=0; i<3; i++){ //new position and Verlet part 2
data[i][n]=dataNew[i][n];
vel[i][n]=vel[i][n] + DeltaT * 0.5*(forceNew[i][n] + force[i][n]) / MASSA;
}
}
totalE=U+K; //total energy
temperature = 2*K/(PARTICELLE*3);
out << t*DeltaT << " " << U  << " " << P  << " " << totalE  << " " << temperature << endl;
}
out.close();
return 0;
}

где моя система находится под потенциалом e ^ -r / r, поэтому у меня есть:

double energy( double r){
return (A*SIGMA*exp(-r/SIGMA)/r);
}
double forza(double r, double h){ //h is for x,y,z
double bubba;
bubba= (A*SIGMA*(exp(-r)*h*(r+1)/(r*r*r)));
return bubba;
}

Спасибо за любую помощь. Я работаю над этим кодом с апреля, и до сих пор у меня нет решения …

редактировать: чтобы быть яснее: КАПИТАЛ условия и DeltaT значения определены в DEFINE