Molecular dynamic, velocity verlet: Kinetic energy divergence

Question

I'm trying to write a simple MD program in C/C++ (I'm used to C but I'm trying to learning C++, so my code is a little "mix"... I know that this is suboptimal and I will move to full C++ as soon as I fully understand it).

Everything seems to run but I have divergences in Kinetic energy, the system does not thermalize and temperature (prop to K) goes from order(10°K) to order(10000°K) in a single step.

I'm working with a low time-step of 0.002 (total time of simulation: 30) so I should not have this enormous error...

This is my code, if something is not clear I can try to explain it better

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;
}

where my system is under a potential e^-r/r, so I have:

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;
}

Thanks for any help. I'm working on this code since April and still I have no solution...

edit: to be clearer: CAPITAL terms and DeltaT are values defined in DEFINE