Kinetic energy goes gradually to zero in my velocity Verlet MD simulation I am trying to figure out what might have caused such a problem: I am simulating a binary Lennard-Jones mixture (two types of atoms interacting via LJ potential) in NVE ensemble with velocity Verlet algorithm. The problem that I have is that, my kinetic energy goes gradually (no sharp change at any single step) to zero, and the potential energy drops gradually (no sharp change at any single step) to a negative value, which is obviously not correct. I printed out the positions and velocities of the atoms, and I find that the velocities go to vanishing values and the positions almost do not change if running very long. The boundary condition is good, so the volume is still fixed. It basically means that in my simulation, the atoms eventually come to a stop. 
It is easy to verify that my velocity Verlet algorithm is correct, but what do you think might be the cause of this kind of weird results? It looks like that the energy is dissipating in my simulation. The code for calculating the potential energy and the force is here, with a cutoff at distance rc, the parameters being the Kob Andersen parameters. Nsmall is the number of smaller atoms, while N is the total number of atoms. For each array (r, f), the first (N-Nsmall) elements are the values for the larger atoms, with the rest Nsmall elements for the smaller atoms.
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;

}
By the way, another question is, should the total energy be strictly constant in a velocity Verlet numerical simulation?
 A: Energy will not necessarily be conserved exactly in any dynamics simulation, because the numerical algorithm is "sampling" the value of some quantity (in your code, apparently the forces fx, fy, fz) during each time step and assuming it remains constant over the whole of the finite length step. In the real world situation, the force is changing continuously during the time step. 
I don't have any experience with molecular dynamics simulations, but this Wikiepdia page has a couple of references to the energy drift problem: https://en.wikipedia.org/wiki/Energy_drift
The Wiki page gives a different (but essentially equivalent) explanation of this behaviour from the first paragraph in this answer: the Verlet algorithm does conserve energy, but the discretized equations that you are solving correspond to a different Hamiltonian from the real-world situation. In general, the difference between the two Hamiltonians will depend on time and the evolution of the numerical solution. Even though it "converges to zero" for each time step considered individually, as the time step size tends to zero, the integral of the error over the whole solution may not converge to zero, or will converge at a slower rate than an error analysis of a single step suggest.
"Conservation of energy" may not be the whole story in any case. For example, there are well known numerical integration methods which conserve energy "exactly" for simple problems like a linear harmonic oscillator, but the constant value of the energy is a function of the time step size, and the incorrect amount of energy in the numerical system results in an error in the oscillation frequency compared with the original differential equation - i.e. the frequency of the numerical solution is not exactly $\omega = \sqrt{k/m}$ but rather $\omega = \sqrt{k/m} + f(h)$ where $h$ is the numerical time step.
