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

Question

Я пытаюсь написать простую 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